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ABSTRACT 

Atmospheric  retrieval  is  the  extraction  of  atmospheric  data  from  spectral  radiance,  as  observed  at 
a  remote  sensor.  In  particular,  consider  the  retrieval  of  temperature  and  humidity  profiles,  and 
aerosol  size  distribution  and  the  scattering  refractive  index  from  long-wave  infrared  and  solar 
radiance  spectra,  respectively.  The  application  of  retrieval,  in  this  report,  primarily  involves 
inversion  of  a  radiative  transfer  equation  (RTE).  However,  due  to  the  ill-posed  nature  of  the 
problem  and  the  inherent  errors  involved,  such  inversions  are  non-trivial.  This  report  presents  a 
combined,  generalised  approach  to  retrieval  via  statistical  inversion,  which  is  derived  in  detail  for 
the  atmospheric  parameters  mentioned  above. 
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Atmospheric  Retrieval  Algorithms  for 
Long- Wave  Infrared  and  Solar  Radiance  Scenarios 


Executive  Summary 


One  of  the  primary  applications  of  electro-optic  sensors  for  the  Australian  Defence  Force  is 
the  accurate  remote  sensing  of  targets  and  background.  Previous  work  of  the  Imaging 
Electro-Optics  Systems  group  at  Defence  Science  and  Technology  Organisation  (lEOS, 
DSTO)  has  shown  that  target  identification  is  sensitive  to  the  influence  and  accurate 
modelling  of  the  intervening  atmosphere.  Consequently,  this  report  focuses  on  acquiring 
an  accurate  representation  of  the  atmosphere  and  on  accounting  for  atmospheric  effects  in 
electro-optic  remote  sensing  scenarios.  The  challenge  is  to  construct  a  valid  description  of 
the  atmosphere  from  observed  spectral  radiance  data. 

The  retrieval  of  atmospheric  parameters,  in  this  report,  is  achieved  through  inversion  of  a 
suitable  approximation  to  the  radiative  transfer  equation  (RTE).  Briefly;  within  the  long¬ 
wave  infrared  spectral  region,  atmospheric  emission  is  the  dominant  contributor  to 
observed  atmospheric  radiance,  characterised  mainly  by  the  vertical  water  vapour  and 
temperature  profiles  of  the  atmosphere.  In  the  solar  spectral  region,  alternatively,  the 
atmospheric  radiance  is  primarily  influenced  by  the  scattering  of  solar  radiation, 
characterised  mainly  by  the  aerosol  and  molecular  properties  of  the  atmosphere.  Eor  these 
two  spectral  regions  then,  appropriate  approximations  to  the  RTE,  dominated  by  either 
atmospheric  emission  or  scattering,  can  be  obtained.  The  constructed  radiative  transfer 
formulae  are  then  inverted  in  order  to  retrieve  the  temperature  and  humidity  profiles  or 
aerosol  size  distribution  and  refractive  index,  for  long-wave  infrared  or  solar  spectral 
sensor  data,  respectively.  Due  to  the  ill-posed  nature  of  the  problem  and  the  inherent 
errors  involved,  such  inversions  are  non-trivial.  The  process  involves  the  application  of 
statistical  regularisation  to  observed  radiance  data  and  a  priori  information. 

This  report  focuses  on  constructing  the  mathematical  and  practical  basis  of  a  generic 
statistical  inversion  method  for  atmospheric  retrieval.  The  purpose  of  developing  a  generic 
base  is  to  demonstrate  how  the  retrieval  procedure  can  be  technically  applied  and  adapted 
to  a  variety  of  scenarios.  In  this  report,  the  process  is  tailored  to  the  scenarios  of  long-wave 
infrared  and  solar  spectral  waveband  sensing.  The  final  result  is  the  construction  of 
appropriate  "retrieval  equations",  obtained  for  the  retrieval  of  temperature  and  humidity 
profiles,  as  well  as  aerosol  size  distribution  and  refractive  index  of  the  atmosphere. 
Extrapolation  of  these  results  is  investigated  with  the  application  to  ground-based  target 
identification,  which  is  found  to  ultimately  rely  on  the  information  available  and  the 
particular  scenario  presented.  Einally,  the  report  concludes  that  further  research  at  DSTO 
would  enable  improved  atmospheric  models  for  Australian  conditions  and  more  accurate 
atmospheric  compensation  in  target  detection  algorithms.  To  this  end,  a  follow-up  report 
will  apply  these  formulae  to  specific  Defence  simulated  and  experimental  scenarios. 
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1.  Introduction 


Atmospheric  retrieval  involves  the  extraction  of  information  about  the  atmosphere  and 
its  constituents  from  remotely  sensed  radiance  data. 

Transfer  of  radiation  through  the  atmosphere  to  a  sensor  is  influenced  by  various 
atmospheric  processes.  The  gases,  aerosols  and  physical  conditions  of  the  atmosphere 
determine  the  extent  to  which  this  radiation  will  be  absorbed,  emitted  and  scattered.  The 
radiance  data  collected  by  a  remote  sensor,  therefore,  contains  information  about  the 
various  processes  and  constituents  of  the  atmosphere.  It  is  ultimately  through 
manipulation  of  this  data  that  atmospheric  characteristics  can  be  retrieved. 

Various  scenarios  and  processes  of  extraction  have  been  studied  in  the  field  of 
atmospheric  retrieval.  This  report  draws  upon  the  plethora  of  research  documentation  in 
the  field  to  present  two  applications  that  have  particular  relevance  to  Defence;  electro¬ 
optic  solar  (ie.  visible  and  short-wave  infrared)  and  long-wave  infrared  sensing.  One  of 
the  primary  applications  of  electro-optic  sensors  in  Defence  is  to  provide  accurate 
remote  sensing  of  targets  and  background.  Previous  work  of  the  Imaging  Electro-Optics 
Systems  group  at  Defence  Science  and  Technology  Organisation  (lEOS,  DSTO)  has 
shown  that  target  identification  is  sensitive  to  the  influence  and  accurate  modelling  of 
the  intervening  atmosphere  [2].  Application  of  atmospheric  retrieval  is  expounded  in 
this  report  in  a  generic  form,  for  these  two  sensing  scenarios. 


1.1  Purpose 

The  purpose  of  this  report  is  to  inform  the  reader  of  the  processes  by  which  to 
technically  perform  atmospheric  retrieval,  via  statistical  inversion.  This  includes  the 
mathematical  construction  of  radiation  transfer  through  the  atmosphere  and  the 
statistical  analysis  necessary  to  retrieve  atmospheric  characteristics  from  observed 
radiance  data.  In  particular,  the  report  will  investigate  two  scenarios;  the  retrieval  of 
atmospheric  temperature  and  humidity  profiles,  and  the  retrieval  of  aerosol  size 
distribution  and  refractive  index.  This  is  made  possible  with  the  data  obtained  through 
passive  sensing  of  long-wave  infrared  and  solar  wavebands,  respectively. 

The  retrieval  processes  for  the  two  different  applications  are  presented  concurrently,  in 
order  to  emphasise  that  atmospheric  retrieval  of  various  phenomena  can  be  conducted 
by  following  a  generic  formula.  Utilisation  of  several  sensors,  of  differing  wavebands, 
could  then  allow  us  to  probe  the  various  aspects  of  the  atmosphere  that  contribute  to  its 
configuration.  Ultimately  this  could  give  us  a  clearer  understanding  of  the  atmosphere 
through  which  we  are  sensing.  One  applications  of  this,  which  is  particularly  relevant  in 
Defence,  is  towards  more  accurate  knowledge  of  the  atmosphere  in  remote  target 
detection. 

The  purpose  of  this  report  is  not  to  validate  a  new  technique  for  atmospheric  retrieval, 
but  rather  to  introduce  an  unfamiliar  reader  to  the  theory  and  mathematics  underlying 
common  atmospheric  retrieval  procedures.  The  point  of  this  is  to  help  readers  apply  the 
mathematics  to  their  own  specific  circumstances.  Application  of  the  formulae  to 
simulated  or  experimental  data  can  be  seen  in  the  referenced  papers  mentioned 
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throughout  the  report.  The  retrieval  formulae  of  this  report  are  well-tested,  applied  and 
documented  in  the  literature,  but  it  is  the  mathematical  development  of  these  equations 
that  is  not  well  documented  in  a  generic  unified  way.  The  importance  of  this  is  two-fold; 
one,  a  report  that  expounds  the  mathematical  derivations  enables  others  to  also 
duplicate  or  follow  the  work  in  this  field;  and  two,  understanding  the  underlying 
assumptions  of  the  resultant  retrieval  formulae  enables  more  accurate  adaptation  of  the 
techniques  to  particular  circumstances  and  scenarios.  Whilst  including  specific  examples 
would  be  beyond  the  scope  of  this  report,  it  is  planned  that  a  less  theoretical  report  will 
follow-up  this  paper  in  the  future,  which  will  apply  these  retrieval  formulae  to  specific 
DSTO  examples. 


1.2  Outline 

This  report  has  be  organised  such  that  the  foundational  mathematics  can  be  accessed  by 
those  that  need  to  adapt  the  forms  to  their  own  needs  (Chapter  2  and  appendices),  but 
also  such  that  the  specific  applications  can  be  examined  without  having  to  delve  into  the 
underlying  mathematics  of  the  statistical  retrieval  (Chapter  3  and  4).  A  summary  of  the 
basic  mathematics  of  the  generic  retrieval  equation  is  presented  in  Chapter  2. 

The  following  chapters  (Chapters  3  and  4)  apply  this  formula  to  the  retrieval  of 
temperature  and  humidity  profiles  from  observed  long-wave  infrared  radiances,  and  the 
aerosol  size  distribution  and  refractive  index  from  observed  solar  radiances.  The 
scenarios  consider  the  radiation  observed  from  ground-based  sensors  directed  with 
zenith  angle  towards  the  sky.  The  sensors  will  detect  radiation  affected  by  various 
processes,  including  radiation  emitted,  scattered  and  absorbed  by  the  atmosphere.  The 
atmospheric  processes  are  mathematically  modelled  according  to  the  application,  and 
then  used  to  retrieve  the  required  information  via  the  retrieval  formula  described  in 
Chapter  2.  More  detailed  explanation  and  derivation  of  the  various  factors  involved  in 
the  retrieval  equations  are  presented  in  the  appendices  (Appendices  A  and  B). 

The  end  of  Chapter  3,  "Application  for  Long-Wave  Infrared  Sensors:  Retrieval  of 
Temperature  and  Humidity  Profiles",  briefly  explores  the  added  application  of  remote 
target  detection.  Within  Defence,  this  is  one  of  the  primary  objectives  of  atmospheric 
retrieval.  As  previously  stated,  knowledge  of  the  atmosphere  and  its  effect  on  the 
radiance,  observed  by  a  remote  electro-optic  sensor,  is  extremely  beneficial  in  accurate 
target  detection.  This  report  will  highlight  the  adjustments  necessary  for  inclusion  of 
target  radiance  in  the  long- wave  infrared  scenario  and  the  added  complications  that  this 
can  produce. 

Finally,  the  report  will  conclude  with  a  summary  of  the  procedures,  an  assessment  of  the 
validity  of  the  applications,  and  suggestions  for  future  progress. 
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2.  Mathematical  Basis  of  Atmospheric  Retrieval 

2.1  Radiative  Transfer  Equation  (RTE) 

The  problem,  as  described  so  far,  is  to  rebieve  values  for  atmospheric  characteristics 
from  experimentally  taken  radiance  measurements.  This  observed  radiance  information 
is  linked  to  the  unknown  characteristics  via  a  Radiative  Transfer  Equation  (RTE). 

The  RTE  describes  the  change  in  spectral  radiance  over  a  distance  {ds)  along  a  path 
through  the  atmosphere. 

Electromagnetic  radiation  transfer  could  be  envisaged  as  an  initial  emission  from  a 
source  (target,  the  Earth  or  the  sun),  which  is  transmitted  through  the  atmosphere  to  a 
sensor.  Various  processes,  involving  atmospheric  gaseous  and  particulate  components, 
affect  the  radiance  along  a  path  from  an  initial  radiator  to  a  sensor.  These  processes 
include  scattering  of  radiation  out  of  the  path,  absorption  of  radiation  along  the  path, 
multiple  scattering  of  indirect  diffuse  radiation  into  the  path,  single  scattering  of  direct 
solar  radiation  into  the  path,  and  emission  of  infrared  radiation  into  the  path.  Thus  the 
RTE,  which  describes  these  processes,  is  expressed  as, 

dL{X)  =  {^)  +  (^l)  +  {X)  -  dL^^^^  {X)  -  {X) ,  eqn(l) 

where  dL{X)  is  the  spectral  change  in  radiance. 

Depending  on  the  applications  of  the  project,  a  RTE  can  be  represented  by  variations  on 
the  above  equation.  As  different  processes  take  dominance  within  different  spectral 
regions  and  circumstances,  the  above  RTE  can  be  simplified  by  removing  less 
consequential  terms. 

In  order  to  appropriately  construct  this  simplified  radiative  transfer  equation,  the  role  of 
the  unknown  characteristics  to  be  retrieved  must  be  considered,  as  will  be  demonstrated 
in  the  application  chapters  to  follow.  This  chapter  will  explore  the  mathematical  basis  of 
the  retrieval  problem  without  going  into  application-specific  detail.  Eirstly,  we  will 
integrate  and  manipulate  the  RTE  into  a  generic,  linearised  form,  as  a  function  of  the 
sensor  radiance.  Then,  using  the  generic  radiance  function  and  other  "information 
equations",  we  will  investigate  the  fundamentals  of  the  statistical  inversion,  as  required 
to  retrieve  atmospheric  characteristics. 


2.2  Discretisation  and  Linearisation  of  the  RTE 

The  radiative  transfer  equation  can  describe  the  transfer  of  electromagnetic  radiation 
through  the  atmosphere  with  an  analytic  function.  It  is  from  this  function,  in  its  various 
forms,  that  we  need  to  extract  information  or  characteristics  of  the  atmosphere. 

In  its  current  form,  the  RTE  (eqn(l))  is  a  function  of  all  atmospheric  processes.  However, 
as  will  be  seen  in  the  following  chapters,  with  specific  application  this  RTE  can  be 
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simplified  to  only  include  the  dominant  processes.  This  is  done  with  the  aim  to  produce 
an  equation  for  the  radiance  at  the  sensor,  as  a  function  of  the  required  atmospheric 
characteristics.  For  the  specific  applications  to  follow  (Chapters  3  and  4),  the  unknown 
atmospheric  characteristics  are  temperature,  humidity,  aerosol  size  distribution  and 
refractive  index.  We  will  represent  these  unknowns  with  the  variable  x . 

Integration  of  the  RTF,  with  respect  to  the  path  length  (ds),  will  give  a  radiance 
equation  for  the  sensor  data  as  a  function  of  the  unknown  atmospheric  characteristics, 

X, 

Aensor  ‘^)  “  dL{^X)  —  (/l,S,x) 

where  .^sensor is  the  observed  sensor  data  at  altitude  z(5),and  L^j^{X,s,x)  is  the 
RTE  integrated  over  the  path  length  to  the  sensor,  as  a  function  of  the  unknown  variable 

X . 

In  these  applications,  as  also  in  other  scenarios,  the  integration  of  the  RTE  is  not 
analytically  possible.  To  solve  this  problem,  we  will  firstly  discretise  the  unknown 
variable  x  into  vector  form,  as  follows. 


2.2.1  Discretisation 

Currently  we  have  assumed  that  x  is  a  continuous  function.  Consider  instead,  that  the 
unknown  variable  is  now  an  unknown  vector,  x . 

For  our  applications,  the  unknown  temperature  and  humidity  profiles  of  the  atmosphere 
are  continuous  functions  of  the  altitude,  T(z)  and  yC>(z) .  Let  us  now  instead  consider  the 
atmosphere  as  'virtually'  constructed  of  discrete  consecutive  layers,  each  of  which  emits 
and  absorbs  electromagnetic  radiation.  We  then  have  layer  altitudes  z.  for  which  we 
can  determine  approximate  temperature  and  humidity  values,  and  hence  a  temperature 
and  humidity  vector,  ze.  x  =  [T^,p^]^,  or  explicitly,  x  =  [7’(Zj),T(z2),...,yC>(Zj),yC>(z2),...]^  . 

Similarly,  the  aerosol  size  distribution  is  a  continuous  function  of  radius,  n{r) .  The 
effects  of  particulate-scattering  processes  on  the  radiation  in  the  atmosphere  are 
fundamentally  determined  by  the  size  of  the  individual  aerosols  and  their  scattering 
cross-section.  Let  us  now  consider  individual  radii  bins  with  which  to  retrieve  a 
discretised  aerosol  size  distribution  vector.  Similarly,  the  refractive  index,  m{A,) ,  can  be 
discretised  into  specific  wavelengths,  corresponding  to  the  radiance  measurements,  ie. 
X  =  or  explicitly,  x  =  . 

Also,  for  the  expedience  of  matrix  formation,  the  observed  radiance  can  be  presented  as 
a  vector  of  wavelengths,  . 

With  appropriate  discretisation  of  the  variables,  and  other  application-specific 
assumptions,  the  RTE  can  be  integrated  along  the  path  to  the  sensor,  and  now  presented 
as  the  radiance  equation. 


L  sensor  ^RTE 


(X). 
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2.2.2  Linearisation 


As  will  be  demonstrated  in  the  applications  to  follow,  the  relationship  between  the 
observed  radiance,  ,  and  the  unknown  atmospheric  characteristics,  x ,  is  often  not 

well  defined  or  simple.  That  is,  the  integrated  RTE  is  not  in  a  form  that  is  practicably 
invertible.  One  solution  is  to  linearise  the  equation  with  a  "first  guess"  of  the 
atmospheric  unknown  variable. 


Let  us  approximate  the  unknown  solution  about  an  initial  estimate,  vector  xO ,  using  a 
multivariate  Taylor's  expansion. 


4ensor(^/)  =  4TE(^.-x)«4TE(2  xO)  +  I  J(2  xO  )  X  -  xO  +  H.O.T.  eqn(2) 


J 


J 


i  V 


J  i 


where 


dL„^p(A.f,x) 

J(A.,xO.)=  J 

J  ^  dx- 


xO 


eqn(3) 


Let  us  assume  for  these  linearised  systems  that  the  higher  order  terms  (H.O.T.)  are 
neglected  without  compromise.  This  will  be  discussed  in  greater  detail  later,  for  the 
individual  applications. 

Hence  we  have  a  discretised  linear  equation  (eqn(2))  that  requires  the  observed  total 
radiance  of  the  atmosphere  data  (2) ,  the  radiance  of  an  atmosphere  with  virtual 

"first  guess"  characteristics  (2,  xO) ,  and  a  calculation  of  how  the  atmosphere  would 

change  with  respect  to  this  virtual  unknown  profile  dL^j^  (2,  xO)/ .  These  parameters 

can  be  calculated  with  the  assistance  of  various  atmospheric  modelling  software 
packages  and  a  priori  information,  as  will  be  discussed  for  the  specific  applications. 


After  discretisation  and  linearisation,  the  general  form  of  the  radiance  equation  as 
expressed  in  equation  (2),  can  be  most  succinctly  expressed  in  a  simple,  generic  matrix 
form, 

L,ensor-LRTE(5^0)  =  J(X-X0) 

eqn(4) 

where  J  is  the  Jacobian  matrix  of  the  derivatives  as  described  in  eqn(3),  and  and 

Lggg(xO)  are  now  vectors  of  specific  wavelengths. 

With  a  generic  equation  (eqn(4))  defined  by  a  simple  matrix  function  of  the  unknown 
atmospheric  vector  and  its  Jacobian  matrix,  it  now  seems  more  feasible  to  retrieve  a 
solution.  Currently,  however,  the  inversion  problem  is  ill-posed;  for  hyperspectral  data 
we  have  an  over-determined  system,  with  more  information  than  unknowns  (and  thus 
no  unique  solution),  and  for  multispectral  data  we  have  an  under-determined  system, 
with  an  inadequate  amount  of  information  for  a  solution. 


5 


DSTO-TR-1879 


It  is  also  evident  that  there  exists  both  measurement  and  estimation  errors.  These 
complicate  the  inversion  procedure  further.  In  order  to  find  the  optimal  solution  with  all 
the  information  available,  whilst  accounting  for  these  errors,  some  form  of  statistical 
regularisation  is  necessary,  tempered  with  knowledge  of  the  physical  reality  of  the 
scenario.  It  is  beneficial  to  firstly  declare  all  information  in  a  generic  form  and  the 
associated  errors,  and  then  extrapolate  to  the  particular  retrieval  scenarios. 


2.3  Information  Equations  for  Retrieval  of  Unknown  Vectors 

The  problem,  as  originally  described,  is  to  retrieve  values  for  atmospheric  unknowns 
from  remote  sensor  radiance  measurements.  This  observed  radiance  data  is  now  linked 
to  the  unknown  atmospheric  vector  via  a  linearised  radiance  equation.  Due  to  the  ill- 
posed  nature  of  the  radiance  equation  and  the  inherent  errors  associated  with  it,  we 
require  more  information  than  that  provided  by  the  sensor.  Statistical  techniques  need  to 
be  used  to  find  an  optimal  solution  that  takes  into  account  all  information  available: 
Consider  these  as  'real'  and  'virtual'  forms  of  information. 

2.3.1  Real  Information 

Real  information  or  data  is  that  which  is  collected  or  measured  from  observed 
phenomena.  Thus  far,  the  sensor  radiances  for  different  wavelengths  and  (later) 

scattering  angles  have  been  considered.  With  any  physical  measurement,  the  equipment 
used  will  always  have  some  degree  of  uncertainty  and  error,  .  Added  to  this, 

the  integrated  radiative  transfer  equation  (eqn(4))  used  to  link  the  observed  real  data  to 
the  unknown  vector  has  inherent  formulation  errors,  due  to  approximations  etc.,  . 

Let  us  attempt  to  account  for  these  errors  in  our  retrieval.  Let, 

=  and  L(x)  =  L^^e(xO)  + J(x-xO)±E„ , 

where  L(x)  is  the  true  value  of  the  radiance;  is  the  observed  data;  is  the 

measurement  {eg.  sensor  instrument)  error;  and  is  the  formula  {eg.  RTE 

approximation  and  linearisation)  error.  So  that, 

Lsensor  =  Lrte (^0)  +  J(x  -  xO)  ±  El ,  whoro  El  =  . 

The  retrieval  procedure,  as  considered  so  far,  has  been  directed  towards  inverting  a 
radiative  transfer  equation  to  obtain  the  desired  atmospheric  variables.  That  is,  the 
primary  information  utilised  has  been  the  observed  total  radiance,  and  the 

radiance  equation  (integrated  RTE)  that  attempts  to  link  the  true  value  of  the  total 
radiance,  L(x),  to  the  unknown  vector  x.  Previously,  and  L(x)  were 

interchangeable,  but  with  addition  of  error  terms,  we  must  redefine  our  equations  as 
above.  The  true  value  of  L(x)  is  unknown  due  to  inherent  errors,  thus  we  are  left  with 
the  question  of  how  to  retrieve  x ,  with  the  observed  radiance  information  ,  within 

the  error  constraints. 
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Let  us  define  the  observed  'real'  information  as  L  * .  The  unknown  variable,  x ,  is 
presently  bound  in  the  integrated  RTE  equation.  Thus,  as  will  be  demonstrated  later,  it  is 
beneficial  to  view  the  current  problem  as  how  to  best  calculate  the  integrated  RTE  from 
the  observed  radiance  data,  in  order  to  extract  x ,  ie. 

L*  =  (xO)  +  J(x  -  xO)  ±  El . 


The  validity  of  the  observed  radiance  information  equation  is,  at  present,  still 
compromised  by  the  assumptions  and  estimations  made  in  the  Taylor's  linearisation. 
The  formula  error  term  in  the  above  equation  can  be  used  to  account  for  the  higher 
order  terms  that  are  neglected  in  this  linearization.  However,  we  have  also  introduced  a 
new  source  of  data,  xO,  without  explanation  of  the  origins  and  inherent  errors 
associated  with  it.  The  "first  guess"  data  xO  is,  in  its  own  right,  a  form  of  information, 
but  since  it  is  derived  from  means  other  than  direct  observation,  let  us  classify  it  as 
'virtual'  information. 


2.3.2  Virtual  Information 

Virtual  information  can  be  obtained  from  statistics  of  the  atmospheric  region, 
extrapolation  of  previous  results,  or  known  constraints  on  the  form  of  the  retrieved  data, 
etc..  Thus  far  we  have  introduced  only  one  form  of  virtual  information;  the  "first  guess" 
profile,  which  is  an  estimation  of  the  retrieved  solution  x , 

x*  =  xO  =  x±Ex. 

where  x  *  is  the  approximated  solution,  often  represented  as  a  first  guess  xO ;  Ex  is  the 
expected  error  in  our  approximation;  and  x  is  the  real  profile  required. 

Linearisation  of  the  integrated  RTE  was  achieved  by  using  a  Taylor's  approximation 
about  a  chosen  parameter  xO.  By  defining  xO  as  the  first  guess  of  x,  and  then 
interpreting  this  first  guess  as  the  real  value  of  x  with  sufficiently  small  error,  we  ensure 
that  the  linearisation  of  the  RTE  is  supportable.  That  is,  with  xO  sufficiently  close  to  x , 
the  linearised  RTE  should  provide  a  reliable  function  about  real  value  x . 


In  the  case  of  the  aerosol  retrieval,  it  is  also  necessary  that  the  retrieved  solution  be  well- 
behaved;  that  is,  smooth  and  continuous.  Let  us  then  introduce  another  virtual 
information  equation  that  will  smooth  the  solution.  As  is  presented  in  Appendix  B.3, 
the  appropriate  smoothing  matrix  B  is  chosen  such  that. 


0*  =  0  =  Bx±Eb. 


2.3.3  Optimisation  of  the  Information  Equations 

All  real  and  virtual  information  about  the  system  can  be  represented  by  a  generic 
function. 
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f*  =  f(x)±Ef . 

eqn(5) 

In  order  to  manipulate  the  various  information  equations,  it  is  recommended  that  the 
functions  f  (x)  be  discretised  and  linearised,  as  already  shown  for  the  observed  radiance 
information  equation  using  Taylor's  approximation  (eqn(4)). 

Assuming  that  xO  is  a  close  approximation  to  x  (with  error  Ex ),  in  generic  terms; 
f(x)  =  f(xO)  + Jj.(x-xO)  +  o((x-xO)^) . 

Therefore  our  "real"  and  "virtual"  information  equations  are; 

Lsensor  =  LrtE  (^0)  +  J(X  -  xO)  ±  El 

xO  =  xO  +  (x  -  xO)  ±  Ex  where  =  dx  /  dx  =  I  and 

0  =  B x0  + J3(x-x0)±Eb  where  Jg=dB/dx  =  B. 


These  equations,  at  first  glance,  may  appear  to  convolute  rather  than  aid  in  the  retrieval 
of  the  unknown  atmospheric  variables  x .  However,  conversion  of  all  information  into 
the  above  forms  is  a  useful  way  to  appreciate  the  amount  and  reliability  of  the 
information  that  we  have  about  x .  From  this  platform,  we  are  able  to  decide  on  a 
method  of  retrieval,  or  more  particularly,  on  a  method  of  statistical  inversion. 


2.4  Statistical  Regularisation  and  Inversion 

We  are  presented  now  with  several  linearised  functions  from  which  to  gain  information 
about  the  unknown  atmospheric  vector,  x  .  How  we  choose  to  use  such  information 
equations  is  dependent  on  the  statistical  technique  best  suited  to  the  application. 

2.4.1  Implementing  the  Inversion 

One  method  of  statistical  inversion  is  to  simultaneously  optimise  all  of  the  information 
equations  with  respect  to  x  .  Each  information  equation  is  considered  optimum  when 
the  error  is  minimised,  such  that  f  *  ^  f  (x)  .  In  terms  of  probability  theory,  we  wish  to 

find  the  x  that  maximises  the  probability  of  finding  {fQ(x),...,f^(x)}  given 
I  fg  f^.  *}  .  This  technique  is  known  as  Maximum  Likelihood  [3] .  The  mathematics  of 

the  statistical  regularisation  technique  employed  here  is  detailed  in  Appendix  A. 

Since,  essentially  we  wish  to  minimise  the  error  variance.  Maximum  Likelihood  is  one 
method  that  facilitates  this.  However,  various  other  statistical  regularisation  techniques 
are  equally  applicable,  and  conclude  with  the  same  final  "retrieval  equation"  forms 
(eqns  (16)  and  (34))  [4]  [5].  The  choice  of  method  depends  on  the  scenario  and 
information  available  [6]  [7]. 


8 


DSTO-TR-1879 


As  described  in  detail  in  Appendix  A,  Maximum  Likelihood  can  be  applied  as  follows: 
By  combining  all  information  equations  (x) ,  such  that  the  errors  are  minimised,  we 
are  able  to  produce  a  formula  for  the  optimum  solution  of  the  atmospheric  variable 
required,  x  (see  Appendix  A.  2); 

=  eqn(6) 

where  are  the  covariance  matrices  of  the  information  equations. 

When  a  specific  application  is  applied,  this  equation  can  be  easily  manipulated  into  a 
retrieval  equation.  That  is,  the  unknown  vector  x  can  be  expressed  as  a  function  of  the 
other  terms  and  variables. 


X  =  fn  (xO),  J,  xO,  Ck ,  b}  . 

eqn(7) 

The  retrieval  of  the  unknown  vector  x  then  becomes  a  search  for  reliable  values  of  the 
terms  in  equation  (7);  the  observed  radiance  the  first  guess  profile  xO,  the 

corresponding  radiance  Lj^^g(xO)  and  its  Jacobian  J ,  and  the  covariance  matrices  C^,. 
This,  in  itself,  is  a  complex  and  error  prone  procedure  with  assumptions,  statistical  data 
and  atmospheric  modelling  packages  as  the  main  tools. 

From  this  phase  onward  it  is  best  to  investigate  atmospheric  retrieval  from  the 
application-orientated  point  of  view.  In  the  following  chapters  we  will  apply  the  above 
mathematical  retrieval  basis  to  the  cases  of  temperature,  humidity,  aerosol  size 
distribution  and  refractive  index  retrieval.  This  will  cover  the  particularities  and 
irregularities  of  the  parameters  needed  and  hopefully  give  a  clearer,  better 
understanding  of  the  abstract  mathematics  presented  above. 


3.  Application  for  "Long-Wave  Infrared"  Sensors: 
Retrieval  of  Temperature  and  Humidity  Profiles 

3.1  Introduction 

This  application  considers  the  radiation  observed  by  a  ground-based  sensor,  sO,  directed 
with  zenith  angle  0  towards  the  sky.  The  remote  sensor  will  detect  radiation  from 
several  sources,  including  that  emitted  from  the  atmosphere,  the  sun  and  any 
obstruction.  All  sources  of  radiation  are  attenuated  as  they  pass  through  the  remaining 
atmosphere,  through  processes  of  absorption  and  scattering.  "The  processes  that  affect 
radiation  along  the  beam  are  scattering  radiation  out  of  the  beam,  absorption  of 
radiation  along  the  beam,  multiple  scattering  of  indirect,  diffuse  radiation  into  the  beam, 
single  scattering  of  direct  solar  radiation  into  the  beam,  and  emission  of  infrared 
radiation  into  the  beam."  [8]. 
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The  constituents  and  physical  condition  of  the  atmosphere  determine  the  significance 
and  magnitude  of  the  radiation  processes  within  specific  spectral  regions.  In  the  visible 
and  short-wave  infrared  range,  radiation  in  the  atmosphere  is  predominantly  effected  by 
atmospheric  absorption  and  the  scattering  of  solar  radiation  into  or  out  of  the  path  of 
observation.  In  the  long-wave  infrared  region,  the  processes  of  atmospheric  absorption 
and  emission,  rather  than  scattering,  are  dominant.  The  ability  of  the  atmospheric 
constituents  to  absorb  and  emit  radiation  is  dependent  on  the  temperature  of  the 
atmosphere  and  the  specific  wavebands  over  which  we  are  sensing.  For  example,  within 
the  8-12gm  region,  water  vapour  (or  the  water  continuum)  is  the  atmospheric 
constituent  that  dominantly  effects  the  radiation.  Thus,  with  the  aim  to  retrieve  both  the 
temperature  and  humidity  profiles  of  the  atmosphere,  one  method  is  to  consider  the  8- 
12pm  long-wave  infrared  radiance  data  from  a  remote  sensor. 

In  the  following  sections  we  will  construct  simplified  radiative  transfer  equations 
specifically  for  the  long-wave  infrared  region.  With  this,  we  will  briefly  revisit  the 
statistical  analysis  of  the  previous  chapter  and  apply  it  directly  to  this  scenario.  Finally, 
we  will  analyse  the  application  to  remote  identification  of  ground-based  targets. 


3.2  Mathematical  Background 

The  first  application  of  atmospheric  retrieval  to  be  considered  involves  solving  for 
unknown  temperature  and  humidity  profiles.  We  will  attempt  to  determine  the  quantity 
of  water  vapour  and  the  effective  temperature  at  any  altitude  in  the  atmosphere. 

The  primary  atmospheric  effects  that  are  influenced  by  temperature  and  humidity  are 
the  absorption  and  emission  of  radiation  by  the  atmospheric  gases.  The  mass  density  of 
water  in  the  air  directly  affects  the  amount  of  electromagnetic  radiation  that  is  absorbed 
by  the  atmosphere.  Temperature,  on  the  other  hand,  is  a  macroscopic  measure  of  the 
random  motion  of  gases  and  particles  that  compose  the  atmosphere. 

In  order  to  retrieve  information  about  the  humidity  and  temperature  profiles,  it  is 
beneficial  to  investigate  the  radiance  measured  in  a  spectral  region  that  is  dominated  by 
the  processes  of  emission  and  absorption  of  radiation  by  atmospheric  gases.  In 
particular,  this  is  served  by  sensing  in  the  long-wave  infrared  region,  8-12gm.  In  this 
scenario,  whilst  the  effects  of  aerosol  scattering  of  radiation  out  of  the  path  must  be 
accounted  for  in  the  transmittance,  the  radiation  observed  from  scattered  solar  radiation 
is  negligible.  (Scattering  of  solar  radiation  into  the  path  contributes  orders  of  magnitude 
less  radiance  than  does  atmospheric  emission,  in  the  8-12pm  band,  according  to 
MODTRAN  execution  [9]).  Thus,  it  is  possible  to  neglect  the  contribution  of  scattered 
radiation  when  constructing  a  radiative  transfer  equation. 

The  atmospheric  constituents  at  any  altitude  will  absorb  and  scatter  out  a  percentage  of 
the  radiation  introduced  to  it;  will  allow  the  remaining  percentage  of  radiation  to  be 
transmitted  onward;  and  also,  will  emit  a  degree  of  thermal  radiation.  The  final  radiance 
that  is  seen  from  a  remote  observer,  looking  either  upwards  or  downwards,  is  the 
accumulated  effect  of  the  transmittance  and  emittance  of  the  atmosphere  on  an  initial 
radiance  {eg.  Earth,  target,  sun). 
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3.2.1  Radiative  Transfer  Equation 

The  radiative  transfer  equation  describes  the  change  in  spectral  radiance  over  the 
distance  ds  along  a  path  (eqn(l)).  For  the  long-wave  infrared  application,  the  RTF  is, 
dL{X)  =  (2)  -  (2)  -  d^^sorp  (^) 

ie.  dL{Z)  =  (2)  -  (2) 

eqn(8) 

where  d4xt  is  the  radiance  absorbed  and  scattered  out  of  the  path 

by  the  atmosphere  (atmospheric  extinction),  and  d4i„ission  i® 

radiance  emitted  into  the  path  by  the  atmosphere. 

CTgxt  (2)  is  the  extinction  coefficient  (measures  the  loss  of  electromagnetic  radiation  due 
to  absorption  and  scattering,  per  unit  distance);  and  B{X)  is  the  blackbody  radiance  of 
the  atmosphere. 

Let  us  define  the  angle  0  as  the  observation  zenith  angle,  dz  =  cos{6)ds ,  where  z  is  the 
altitude;  and  let  us  define  the  incremental  optical  depth  as  dr^^^  =-co?,{6)o'^^^ds  . 

For  down-welling  radiation,  optical  depth  is  defined  as, 

^ext  is)  =  cos(6')£  (2, 5  ')ds ' . 

Thus,  the  RTF  (eqn(8))  becomes, 

cos(6')  cos(6') 

In  order  to  describe  the  observed  radiance  data  as  a  function  of  the  atmospheric 
processes,  let  us  integrate  the  RTF  above.  For  a  ground-based  sensor,  integrating  over 
the  entire  path  length  from  space  to  sO  gives  the  radiance  equation, 

^  -r' 

COS(C7j  tot^  ' 

which  is  simply  =  4ni,ial^''««'^atmosphere  +  Atmosphere  '  eqn(10) 

The  transmittance,  trans{X),  is  defined  as  the  fraction  of  radiation  allowed  to  pass 
through  a  portion  of  the  atmosphere,  without  being  absorbed  or  scattered  out.  For  the 
downwelling-radiation  case  {eg.  ground  sensor  to  space). 
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trans{X,s) 


L{X,s) 

Anitial("^’°°) 


L 


initial 


r(/l,.y) 


initial 


exp 


eqn(ll) 


If  we  were  to  position  the  sensor's  line  of  sight  (path)  directly  at  the  sun,  then  the  initial 
radiation,  (A,) ,  would  be  the  initial  solar  radiation.  Let  us  assume  this  is  not  the 
case  in  our  scenario,  so  that  there  is  no  initial  radiation  and  the  sensor  records  only 
atmospheric  effects.  Thus, 


-z~' 

(A,  sO)  =  4„  (1)  =  S(l)  e'”'" Vr  ■  =  dtrans . 

eqn  (12) 

The  integrated  RTE,  radiance  equation  (eqn(12)),  is  now  simplified  as  a  function  of 
blackbody  atmospheric  radiance  and  transmittance. 

The  analytical  form  of  Planck's  blackbody  equation  B(T(s),  A,) ,  is  well  known  and 
defined  as  a  function  of  temperature  and  wavelength  alone,  and  thus  is  a  function  of  the 
altitude  in  so  much  as  the  temperature  varies  through  the  atmosphere. 

The  transmittance,  trans(/i,s) ,  is  not  a  well-behaved  or  differentiable  function.  The 
extinction  coefficient  a^^^(A,,s) ,  from  which  it  is  derived,  is  dependent  on  the 
composition  of  the  atmosphere  at  a  particular  height,  and  its  corresponding  extinction 
properties.  The  extinction  coefficient  and  transmittance  are  directly  dependent  on  the 
water  composition  or  absolute  humidity  of  the  atmosphere  at  any  altitude.  The 
temperature  is  indirectly  related  to  the  transmittance,  as  a  mutual  consequence  of  the 
altitudinal  composition.  In  an  unfamiliar  atmosphere  the  exact  composition  and 
properties  of  all  elements  in  the  atmosphere  is  unknown.  The  long-wave  infrared 
waveband  used  for  this  retrieval  is  dominantly  influenced  by  water  vapour  content  and 
temperature,  which  enables  us  to  use  model  approximations  of  the  other  constituents 
without  too  great  an  error  or  complications.  Thus  transmittance  can  be  described  by  the 
temperature  and  humidity  profiles  of  the  atmosphere.  However,  since  the  extinction 
coefficients  are  ultimately  only  known  through  empirical  observation  (data  tables),  we 
do  not  have  transmittance  as  a  well-behaved,  continuous  function  of  temperature  and 
humidity. 

In  summary,  the  radiance  equation  above  is  not  presently  in  an  invertible  form. 
Transmittance  is  an  ill-posed  function  of  temperature  and  humidity,  as  well  as  various 
other  factors  that  are  unknown  for  the  atmosphere  of  interest.  We  can  approximate  the 
other  characteristics  to  an  acceptable  degree,  but  require  further  approximations  before 
inversion  retrieval  can  be  performed. 


3.2.2  Assuming  a  Layered  Atmosphere 

Consider  the  unknown  atmosphere  as  structured  with  consecutive,  homogeneous 
altitudinal  layers.  Each  layer  contributes  to  the  extinction  and  emittance  of  radiation  in 
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the  path  to  the  sensor.  It  is  through  introduction  of  atmospheric  layers  (and 
linearisation)  that  we  will  solve  the  problem  of  ill-posed  transmittance  functions,  and 
thus  the  retrieval  of  temperature  and  humidity  for  specific  altitudes. 

-t' 

4,,  {X,sO)  =  f  B(A,  T{t  '))  dr ' 

cos(6») 

Extraction  of  the  unknown  temperature  and  humidity  profiles  from  the  integral 
equation  above  is  analytically  unachievable.  We  do  not  have  analytic  equations  for  the 
transmittance  or  optical  depth  as  functions  of  the  unknown  profiles.  Empirical  data 
tables,  however,  exist  for  the  extinction  coefficients  under  specific  conditions.  Therefore, 
if  we  consider  the  atmosphere  to  be  homogeneous  over  chosen  altitudinal  intervals,  we 
can  establish  values  of  optical  depth  for  specific  values  of  temperature,  humidity  and 
altitude,  t{T.  ,  p.  ,X,s)  .  We  consider  optical  depth  rather  than  transmittance  here,  since  it 
is  proportional  to  altitude,  and  thus  more  accurately  discretised  with  respect  to  altitude. 

Therefore,  let  us  firstly  define  the  unknown  temperature  and  humidity  functions  as 
atmospheric  profiles  in  vector  x , 

where  T  =  [7’(z,),7’(z,),...f  and  p  =  [p(Z|),p(z, 


and  let  us  assume  that  the  optical  depth  can  be  described  by  the  temperature  and  water 
content  of  the  atmosphere  alone,  t{Z, z.)  =  t(X, z. , T, p)  . 


So  that  we  can  now  describe  the  radiance  equation  as  a  discrete  summation  over 
atmospheric  layers,  as  follows, 

-r(T.,p.) 


-1  COS(^) 


,=1  cos(6>) 


Ar(T,,p,), 


or  simply. 


Arte  sO)  =  ^  z. ,  T,  p)  trans{A,  z. ,  T,  p) . 


i=i 


eqn(13) 


eqn  (14) 


The  radiance  equation  (eqn(14))  now  presents  the  sensor  observed  radiance  as 
constructed  from  the  sum  of  the  emitted  radiation  from  each  layer  in  the  atmosphere 
(which  are  in  turn  transmitted  through  the  rest  of  the  atmosphere  to  the  sensor). 


Each  atmospheric  layer  is  composed  of  varying  percentages  of  gases  and  particles, 
including  water  vapour.  It  is  the  thermal  properties  and  condition  of  these  constituents 
that  determine  the  temperature  at  any  altitude.  The  discretised  radiance  equation 
(eqn(14))  assumes  homogeneity  within  the  constructed  layers  of  the  atmosphere.  Eor  our 
calculations,  therefore,  there  is  an  assumption  that  the  temperature  (and  thus  Planck's 
radiance)  and  the  constituent  concentrations  are  constant  within  each  layer.  This 
discretisation,  then,  is  only  valid  for  sufficiently  small  layers.  Consideration  as  to  the 
size  of  the  layer  must  include  the  atmospheric  conditions  and  the  positioning.  Eor 
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example,  layers  closer  to  the  ground  will  need  to  be  narrower  than  those  at  higher 
altitudes,  due  to  the  relatively  steep  change  in  temperature  with  altitude  near  the  Earth's 
surface. 


3.2.3  "First  guess"  Temperature  and  Humidity  Profiles 

The  blackbody  radiance  emitted  by  each  atmospheric  layer  is  a  well-posed  function  of 
r(z.)  and  yC>(z.),  however,  the  transmittance  trans{A,z.)  and  optical  depth  r(/l,z.)  are 
not.  An  analytic  statistical  inversion  of  the  current  radiance  equation  (eqn(14)),  then,  is 
not  possible. 

One  solution  is  to  linearise  the  equation  with  a  "first  guess"  of  the  atmospheric 
conditions.  In  this  way,  the  transmittance  for  the  "first  guess"  atmosphere  can  be 
empirically  sought,  and  our  equation  becomes  an  explicit,  invertible  function  of  the 
unknown  temperature  and  humidity  vector. 


Let  us  approximate  the  radiance  equation  about  an  initial  profile,  xO  =  [TO^,  pO^]^, 
using  a  multivariate  Taylor's  expansion; 


4ensor(^y)~4TE(^y.xO)  +  X 


^  (  dL 


n=l 


a) 


RTE 


5T„ 


(T„-T0„) 
xO  y 


P=1 


(Pp-pOp) 

xO  y 


where  T^ensorC"^)  is  the  observed  data,  and  TR.pg(/l,xO)  is  the  radiance  equation  (eqn(14)) 
as  a  function  of  "first  guess"  temperature  and  humidity  xO . 

In  matrix  form,  the  above  linearised  equation  becomes, 

Lse„sor-LRTE(XO)«  J(X-XO) 


where  J  is  the  Jacobian  of  the  radiance  equation  (eqn(14)). 


5Trte(^,) 


J,  I 


dX; 


eqn(15) 


xO 


and  i  =  n  +  p  is  the  number  of  the  atmospheric  layer  being  described,  and  j  is  the 
number  of  the  radiance  wavelength. 

Hence  we  have  a  discretised  linear  equation  that  requires  the  observed  total  radiance, 
^sensor '  t^e  radiance  values  of  the  atmosphere  at  varying  altitudes  for  selected  "first 

guess"  atmospheric  conditions  of  temperature  and  humidity.  Now  that  we  have  an 
expression  in  the  generic  form  of  equation  (4),  as  required  in  Chapter  2,  we  can 
investigate  the  inclusion  of  error  and  the  other  information  equations  and  from  thence, 
the  statistical  inversion. 


+H.O.T. 
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3.2.4  Temperature  and  Humidity  Profile  Retrieval  Equation 

For  the  case  in  which  we  are  required  to  retrieve  temperature  and  humidity  profiles, 
smoothing  (Appendix  B.3)  is  not  necessary.  Thus  we  only  include  two  information 
equations  in  our  minimisation;  the  linearised  radiance  equation  above,  and  the  "first 
guess"  approximation  of  the  profiles,  xO .  By  combining  the  information  equations  so 
that  the  errors  are  minimised,  we  are  able  to  produce  a  formula  for  the  optimum 
solution  to  the  atmospheric  parameter  required,  x .  Appendix  A.2.1  has  details  of  the 
minimisation  technique  used  and  the  specific  application  to  temperature  and  humidity 
profile  retrieval. 

By  substituting  fQ(x)  =  Lj^^g(xO)  and  fj(x)  =  x  into  the  minimisation  equation  in 

Chapter  2  (eqn(6)),  as  seen  in  Appendix  A.2.1,  we  can  rearrange  and  manipulate  our 
minimisation  equation  into  the  form  below. 

This  is  now  the  referral  equation  by  which  we  can  calculate  retrieved  values  for 
temperature  and  humidity. 


eqn(16) 

where  x  is  the  unknown  temperature  and  humidity  profile;  is  the  measured 

radiance  from  the  sensor;  xO  is  the  approximate  estimate;  Lj^^g(xO)  is  the  radiance  of 

this  estimate;  J  is  the  Jacobian  (eqn(15));  and  the  C  are  the  covariance  matrices  of  the 
errors,  as  described  in  Appendix  B.2. 


This  completes  the  mathematical  background  of  the  atmospheric  retrieval  in  the  long¬ 
wave  infrared  region.  We  have  developed  a  specific  retrieval  equation  that  can  be  used 
to  solve  for  the  temperature  and  humidity  profiles  of  an  unknown  atmosphere  (eqn 
(16)).  From  this  stage  onward,  the  specific  application  and  scenario  must  be  considered, 
and  the  use  of  comprehensive  modelling  packages  becomes  prominent.  The  following 
section  will  summarise  the  findings  above  and  give  a  brief  overview  of  the  practicalities 
involved  in  its  application. 


3.3  Temperature  and  Humidity  Retrieval:  Summary  and  Practicalities 

In  summary,  we  consider  the  scenario  in  which  we  are  ground  sensing  the  radiance,  at 
long-wave  infrared  wavelengths,  which  has  been  emitted  and  transmitted  by  the 
atmosphere.  We  are  interested  in  retrieving  values  for  the  temperature  and  humidity  of 
the  atmosphere  from  this  observed  radiance  data. 

Due  to  a  lack  of  a  well-defined  transmittance  function,  as  previously  explained,  the 
spectral  data  from  a  remote  sensor  must  be  supplemented  by  a  reasonable  "first  guess" 
of  the  temperature  and  humidity,  from  which  we  can  extrapolate  the  true  solution.  An 
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understanding  of  the  errors  involved  in  the  measurement  and  the  retrieval  procedures  is 
also  necessary  for  accurate  atmospheric  retrieval. 

Through  statistical  regularisation  of  the  information  available,  the  following  equation 
can  be  created,  via  which  we  can  retrieve  the  temperature  and  humidity  of  the 
atmosphere  at  consecutive  altitudes. 


x=xo+(rc,|i+qi)‘rc,i  (l 


r' 


,(x0)) 


where, 

X  is  the  vector  of  unknown  temperature  and  humidity  variables  to  be  retrieved, 

X  =  [r(zj ),  r(z2 ),...,  r  (z^ ),  p{z, ),  p{z^ ),...,  p{zp)f ; 

L^ensor  ^^e  vector  of  measured  radiance  values  from  the  sensor, 

^sensor  t "^sensor  X -^sensor  ("^2  X  ***  5 -^sensor  ' 

xO  is  the  vector  of  first  guess  approximations  of  the  unknown  x , 
xO  =  [TO(Zj ),  T0{z^ ),...,  T0{z^  ),  /?d(zj ),  pO{z^), pO{Zp  )]^ ; 

L^2^g(x0)  is  the  vector  of  radiance  created  by  simulating  an  atmosphere  with  xO  profile, 
using  eqn(14), 

Lrte  (xO)  —  [^RXE (A  ’  xO),  (^2 ,  xO), . . . ,  ,  xO)]  , 


J  is  the  Jacobian  matrix  of  the  radiance  equation,  eqn(15).  The  derivatives  are  calculated 
with  respect  to  the  unknown  variables  and  then  calculated  for  xO , 


54te(^,) 


5T„ 


and 


54te(^/) 


'j,p+N 


xO 


xO 


for  J  =  1...M  wavelengths,  n  =  1...N  temperature  altitudinal  layers,  and  p  =  \...P 
pressure  altitudinal  layers; 


and  C  are  the  covariance  matrices  that  account  for  inherent  errors  in  the  process. 
Cm  and 

Cbx  =i[c,,(r(z,)),...,c,,(r(z^)),c,,(p(z,)),...,c,,(p(z,))f . 


We  initially  assume  that  the  covariance  matrices  are  diagonal  (that  is,  errors  are 
independent).  However,  iteration  during  the  retrieval  process  can  allow  for  variation,  as 
will  be  explained  in  section  3. 3.1. 2. 


Establishing  the  "optimal"  inversion  equation  (eqn(16))  is  the  first  step  in  the  retrieval 
procedure.  Following  this,  the  individual  components  must  be  correctly  calculated  and 
the  incorporated  errors  must  be  estimated  to  an  acceptable  degree.  This  involves  the 
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composition  of  non-trivial  matrices  and  many  strictures  imposed  on  the  processes 
involved,  as  will  be  discussed  below.  (More  detailed  descriptions  of  the  Jacobian  and 
Covariance  matrices  are  presented  in  Appendices  B.l  and  B.2.) 


3.3.1  Practicalities 

3. 3.1.1  Atmospheric  Models 

For  the  atmosphere  under  observation,  the  exact  gas  and  aerosol  composition  is 
unknown.  Atmospheric  modelling  packages,  like  MODTRAN,  offer  several  models  that 
represent  various  'typical'  atmospheres.  These  model  atmospheres  are  composed  of 
specific  combinations  and  quantities  of  water,  carbon  dioxide,  ozone,  pollution,  etc.  and 
refer  to  average  climatological  temperature  profiles.  The  optical  properties  of  the 
atmospheric  constituents  have  been  empirically  tested  for  each  model  in  MODTRAN 
and  the  extinction  coefficients  are  data-based  within  the  program. 

The  "first  guess"  temperature  and  humidity  profiles,  xO,  can  be  calculated  through 
atmospheric  simulation  packages,  such  as  MODTRAN  and  FASCODE  [10],  or  taken 
from  meteorological  data  available  for  the  region.  Depending  on  the  sensing  scenario, 
an  appropriate  model  can  be  chosen  to  represent  the  atmosphere.  The  temperature  and 
humidity  profiles  of  this  model  atmosphere  can  then  be  used  as  the  "first  guess"  profile 
xO .  The  linearisation  of  the  radiance  equation  that  occurs  in  the  retrieval  calculations 
requires  this  first  guess  to  be  adequately  close  to  the  true  value  of  x  .  This  can  lead  to 
difficulties  if  the  location  that  is  being  observed  is  foreign  and  unmodelled.  One  method 
of  calculation  of  initial  guess  profiles  from  ground  based  sensors,  as  developed  by 
Theriault  et  al.  [11]  [12],  involves  manipulation  of  radiance  data  collected  at  highly 
absorbent  frequencies  eg.  for  CO2.  Alternatively,  statistical  or  climatological  profiles  can 
be  used  as  an  appropriate  "first  guess"  [13]. 

From  the  initial  approximation  of  the  unknown,  xO ,  the  values  for  the  virtual  radiance 
L^.j.g(x0)  and  the  Jacobian  matrix  J  must  be  calculated.  Unlike  the  observed  radiance, 
the  radiance  that  corresponds  to  the  first  guess  profile  cannot  be  measured  directly.  A 
suitable  RTE  (eqn(14))  must  be  used  to  calculate  Lj^^g(xO)  from  the  'virtual' 

atmospheric  profile  xO .  Options  for  the  user  to  input  their  own  temperature  and  gas 
profiles  are  available  in  MODTRAN  or  EASCODE,  and  it  is  through  this  that  values  for 
the  virtual  radiance,  L^^g(xO) ,  can  be  created.  The  inclusion  of  aerosol  extinction  and 
the  effect  of  other  particulate  and  gaseous  constituents  in  the  accurate  formation  of  the 
RTE  and  retrieval  work  are  often  heavily  based  on  a  priori  model  atmospheres  such  as 
those  in  MODTRAN. 

Similarly,  the  Jacobian  of  the  radiance  equation,  J ,  requires  a  simulation  package  (or 
data-base)  to  model  the  coefficients  of  an  unknown  atmosphere.  As  previously 
explained,  we  cannot  form  an  analytic  function  of  the  transmittance  with  respect  to 
temperature  and  humidity.  Thus  direct  calculation  of  the  derivatives  of  the  radiance 
equation  is  impossible.  The  Jacobian  matrix,  therefore,  is  constructed  by  approximating 
the  derivatives  by  assessing  the  effect  of  small  linear  changes  in  the  temperature  to  the 
observed  radiance.  The  radiance  is  calculated  for  the  "first  guess"  profile  and  the 
subdominant  constituents  are  approximated  by  the  appropriate  atmospheric  model,  in 
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order  to  create  estimates  of  the  transmittance  as  required  in  the  formula.  More  detailed 
explanation  of  this  process  is  provided  in  Appendix  B.l. 


3. 3.1.2  Newtonian  Iteration  of  the  Solution 

The  form  of  our  retrieval  equation  (eqn(16))  is  such  that  an  iterative  process  can  be 
applied  to  acquire  more  accurate  results.  We  can  substitute  the  resultant  x  back  into  the 
equation  as  the  new  xO  value  and  recalculate.  This  practice,  however,  may  lead  to 

misleading  results  if  the  Jacobian  J  and  first  guess  radiance  Lj^.j,g(xO)  are  not  also 
updated  with  the  new  xO .  The  majority  of  the  convergence  towards  the  real  solution 
occurs  within  the  first  iteration.  The  covariance  matrix  ,  also,  can  be  iteratively 
calculated  with  statistical  regularisation.  In  this  way,  a  better  adjustment  for  error, 
which  accounts  for  possible  correlations,  can  be  acquired.  For  further  information  on 
this  approach  to  the  covariance  matrices  see  Rodgers  [14]. 


3. 3. 1.3  Sources  of  Error 

The  accuracy  of  L^.pg(xO)  relies  on  the  suitability  of  the  modelling  software  program 
used  and  the  accuracy  of  the  "first  guess".  The  errors  inherent  in  the  calculation  of 
L^.j.g(xO)  therefore  come  from  two  sources;  "first  guess"  error  variance  of  the 
temperature  and  humidity  profiles,  and  error  in  the  model  representation 
approximations  of  all  other  atmospheric  contributions.  The  first  of  these  errors  is 

accounted  for  in  the  "first  guess"  error  covariance  matrix  of  the  retrieval  equation,  . 
This  allows  the  solution  to  vary  within  defined  limits.  The  error  in  the  atmospheric 
model  RTF  is  hopefully  diminished  by  the  appropriate  choice  of  spectral  range.  The 
dominant  radiative  processes  in  the  long- wave  infrared  range  should  be  summarised  by 
the  temperature  and  humidity  profile  of  the  atmosphere.  The  scenario  and  specific 
circumstances  of  the  application  will  define  whether  this  is  accurate  enough  for  the 
project. 

The  Jacobian  matrix,  J,  is  limited  by  the  errors  above,  and  also  by  the  practical 
approximation  that  is  required  for  its  calculation.  The  analytic  derivatives  of  the  RTE 
cannot  be  calculated  due  to  the  discrete  nature  of  the  transmittance.  Thus  the  "small 
linear  change"  of  the  radiance  equation  with  respect  to  x ,  as  explained  in  Appendix  B.l, 
will  unavoidably  produce  error.  This  is  one  source  of  error  that  must  be  considered 

when  constructing  the  radiance  error  covariance  matrix  of  the  retrieval  equation,  Cg, . 

This  matrix  attempts  to  account  for  the  errors  in  the  first  RTE  information  equation, 
which  includes  the  error  involved  in  the  linearisation  of  the  RTE  and  in  any  of  the  terms. 
The  radiance  error  covariance  matrix  also  must  account  for  the  error  in  observational 
measurement.  This  should  be  the  known  variance  error  of  the  sensor. 
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3.4  Atmospheric  Retrieval  for  Target  Detection 

One  motivation  for  the  retrieval  of  atmospheric  characteristics  is  the  need  for  better 
representation  of  atmospheric  influence  in  remote,  passive  target  detection.  We  will  now 
investigate  the  scenario  of  a  remote  hyperspectral  sensor  collecting  radiance  data,  in  the 
long- wave  infrared  region,  from  a  target  on  the  ground.  The  result  of  this  is  to  hopefully 
determine  the  characteristics  and  identify  this  target,  through  consequent  calculation  of 
its  emissivity  and  temperature. 

Previous  study  by  the  lEOS  group  at  DSTO  Edinburgh  (South  Australia),  has  indicated 
that  the  accurate  detection  of  remote  targets,  in  this  scenario,  requires  adequate 
representation  of  the  atmospheric  characteristics  [13].  It  was  found  that  even  1% 
improvement  of  the  temperature  profiles  and  5%  improvement  of  the  humidity  profile 
gained  substantially  better  results  for  retrieval  of  temperature  and  emissivity  of  a  target. 
It  is  with  this  motivation  that  achieving  more  accurate  temperature  and  humidity 
profiles  of  the  atmosphere  between  the  target  and  sensor  is  considered  beneficial. 

The  scenario  presented  here  differs  from  the  previous  section  in  that  we  are  now 
interested  in  the  upwelling  radiation  observed  by  a  remote  sensor.  The  radiative  transfer 
equation  must  accommodate  the  inclusion  of  an  initial  radiance  and  adjust  for  this  new 
orientation. 


3.4.1  Modelling  Upwelling  Radiative  Transfer 

The  application  described  in  this  chapter,  "Application  to  Retrieval  of  Temperature  and 
Humidity  Profiles",  is  constructed  with  the  aim  towards  characterising  the  atmosphere 
without  reference  to  an  initial  emitter.  The  scenario  was  described  as  a  ground-based 
sensor  directed  into  the  sky.  Let  us  now  re-adjust  our  perspective  so  that  we  have  a  sky- 
based  sensor  directed  at  the  Earth's  surface.  The  first  problem  that  must  be  dealt  with, 
then,  is  to  investigate  the  remotely  observed  radiance  over  ground  surface  that  is 
known.  The  method  for  doing  this  is  dependent  on  the  specific  scenario  of  the 
experiment. 


In  section  3.2.1  we  presented  the  observed  sensor  radiance  as  a  function  of  the  initial 
radiation,  atmospheric  radiation  and  atmospheric  transmittance. 


“^sensor  "^initial^^^^'^atmosphere  "^atmosphere  ' 


Let  us  present  this  now  in  the  appropriate  orientation-  for  upwelling  radiation  as 
observed  by  a  remote  sky-based  sensor  at  height  z(S) ; 

Censor  i^,S)  =  4^^,  {A,  sO)  trUHsl^  (A)  +  4m  W  ■ 


”1" 

where  4tm  up-welling  atmospheric  radiation  and  transmittance. 

The  path  of  the  radiation  is  now  between  the  sensor  S  and  the  ground  sO ,  rather  than 

space.  Total  transmittance  is  now  defined  as,  trans^(A,  S)  =  exp  |-J  ^  {A,  s  ')ds  '| . 
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3A.1.1  Ground  Surface  Radiance 

In  the  previous  sections  the  initial  radiation  term  was  not  applicable.  We  now  consider 
the  case  in  which  the  initial  radiation  of  the  Earth  is  observed  and  thus  included  in  our 
radiance  equation.  The  initial  Earth  radiance  is  a  result  of  the  emitted  thermal  radiation 
of  the  Earth  (and  target)  and  the  reflected  solar  and  atmospheric  radiation  from  the 
Earth's  surface. 


"^earth  "^emitted  "^reflected 

4-  =  p„,.  (X)(eI^(X,sO)  +  eIjX,sO)).  eqnd  7) 

The  emissivity,  s(/l) ,  is  a  measure  of  how  effectively  a  material  can  emit  radiation; 
defined  as  the  ratio  of  emitted  radiation  from  the  source  (Earth)  to  the  emitted  radiation 
from  a  blackbody  of  the  same  temperature.  It  is  specific  to  the  material  investigated  and 
a  function  of  wavelength.  Blackbody  radiation  B{A,T)  is  solely  a  function  of 
wavelength  and  temperature.  The  reflectivity,  p(A,) ,  of  a  solid  opaque  object  in  thermal 

equilibrium  with  its  surroundings  is  (l-f(/l)).  Lloi^r/atmi-^AO)  is  the  downwelling 

radiance,  contributed  by  the  sun  or  the  atmosphere,  as  measured  at  the  Earth's  surface, 
sO. 

In  order  to  retrieve  information  about  the  atmosphere  from  the  upwelling  radiation 
measured  at  the  sensor,  we  firstly  need  to  know  the  emissivity  and  temperature  of  the 
ground  surface.  This  leaves  us  with  a  circular  problem;  the  purpose  of  target  detection 
analysis  is  to  retrieve  (ground)  target  emissivity  and  temperature,  which  are  themselves 
necessary  for  the  retrieval  equation. 

Atmospheric  retrieval  is  a  useful  tool  in  target  detection  for  acquiring  a  better 
approximation  of  the  atmosphere  than  a  "first  guess"  simulated  model.  Let  us  assume, 
in  the  following  section,  that  there  is  a  region  of  ground  surface  sufficiently  close  to  the 
target,  which  has  a  recognisable  emissivity  signature.  By  considering  the  atmospheric 
characteristics  retrieved  nearby  the  target,  this  will  hopefully  give  a  better 
approximation  of  the  atmosphere  in  the  line  of  sight  to  the  target. 


3.4.2  Retrieval  of  "Surrounding"  Atmospheric  Characteristics 

Let  us  attempt  to  retrieve  the  atmospheric  characteristics  by  analysing  the  radiance 
observed  in  the  area  surrounding  the  target,  for  which  the  surface  material  is  familiar. 
The  simplest  case  of  this  would  be  if  the  ground  surface  were  approximately  a 
blackbody  ( ^  =  1 )  or  greybody,  for  example,  the  ocean,  which  has  almost  constant 
emissivity  [14].  Other  scenarios  include  familiar  ground  surfaces  for  which  the 
emissivity  is  known,  for  example  concrete  or  desert.  In  order  for  the  assumption  of  a 
priori  information  about  the  emissivity,  it  may  be  necessary  for  all  the  suitable  data  from 
the  surrounding  area  to  be  collated  and  used  as  an  average  radiance  measure  of  the 
average  surrounding  atmosphere  and  surface  properties. 
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The  average  temperature  of  the  surrounding  ground  surface,  T (sO) ,  is  then  included  as 
an  element  of  the  unknown  vector,  x ,  with  extra  elements  being  added  to  the  first  guess 
vector,  xO ,  and  Jacobian,  J .  This  is  not  done  with  the  final  aim  to  find  an  accurate 
temperature  for  the  average  surface  T (sO) ,  but  to  ensure  the  use  of  an  accurate  ground 
temperature  in  the  retrieval  of  the  average  atmospheric  properties  of  the  surrounding 
area.  This  scenario  requires  Earth  radiance  terms,  with  known  emissivity,  to  be  included 

in  the  relevant  RTEs  that  are  needed  to  calculate  L^.j.g(xO)  and  J.  With  this 
information,  the  retrieval  equation  below,  unchanged  from  the  previous  section,  can  be 
used  to  retrieve  the  temperature  and  humidity  profiles  for  the  atmosphere  surrounding 
the  target; 


x=xo+(rc,|i+q;)  ‘rc,i  (l 


r' 


,(x0)) 


where  now 

eqn(18) 

for  4tmW= 

j=sO 

and  =  trans^{X,T,^),  (see  section 4.2.2.1). 


The  atmospheric  properties  of  temperature  and  humidity  retrieved  here  can  then  be 
used  to  determine  approximate  transmittance  and  atmospheric  radiance  values  for  the 
line  of  sight  of  the  sensor  to  target.  These  can  then  be  used  to  calculate  the  emissivity 
and  temperature  of  the  ground  target. 

This  method  will  invariably  add  more  uncertainty  to  the  retrieval  of  the  atmospheric 
temperature  and  humidity.  The  ground  surface  radiance  is  the  dominant  contributor  in 
the  long-wave  infrared  range  that  we  are  examining.  The  inclusion  of  correspondingly 
large  Jacobian  terms,  due  to  this,  will  cause  less  sensitive  retrieval  of  atmospheric 
profiles.  Thus  this  approach  cannot  retrieve  the  atmospheric  characteristics  with  as 
much  accuracy  as  the  case  where  the  surrounding  ground  temperature  is  known. 


3.4.3  Alternative  Technique 

Accurate  modelling  of  the  atmosphere  is  essential  for  good  estimation  and  retrieval  of 
ground  characteristics,  and  thus  target  identification.  As  such,  any  improvement  on  the 
assumed  or  "first  guess"  atmosphere  model  is  worth  investigating.  The  accuracy  of 
retrieval  via  statistical  inversion  is  dependent  on  the  information  available  and  the 
scenario  involved,  as  is  presented  above. 

In  regards  to  target  detection  for  unknown  or  foreign  Earth  surfaces,  the  inversion 
method  can  facilitate  the  retrieval  of  atmospheric  characteristics  for  only  limited 
scenarios.  An  alternative  approach  is  to  use  a  numerical  optimisation  method.  This 
involves  repeatedly  forward-modelling  the  scenario  with  different  atmospheres  until  the 
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"best"  model  solution  is  found.  An  example  of  this  technique  is  given  in  [17]  and  [18]. 
This  approach  is  best  suited  to  the  classification  of  targets,  rather  than  accurate  retrieval 
of  atmospheric  characteristics.  Thus  the  choice  of  retrieval  method  is  dependent  on  the 
final  aim  of  the  project. 


4.  Application  for  "Solar"  Sensors:  Retrieval  of 
Aerosol  Size  Distribution  and  Refractive  Index 

4.1  Introduction 

In  the  long- wave  infrared  region,  the  temperature  and  water  content  of  the  atmosphere 
can  be  retrieved  due  to  the  dominant  influence  of  absorption  and  emission  on  the 
radiation  in  an  observation  path.  In  the  solar  spectral  region,  however,  the  effects  of 
atmospheric  thermal  emission  are  surpassed  by  the  influence  of  scattering  of  solar 
radiation  into  and  out  of  the  path,  by  the  atmospheric  constituents.  Aerosol  size  and 
refractive  index  influence  the  degree  of  scattering  of  solar  radiation  by  an  atmospheric 
particle.  Thus  through  manipulation  of  the  observed  atmospheric  radiance  in  the  visible 
and  short-wave  infrared  spectral  regions,  we  can  potentially  gain  information  about  the 
atmospheric  aerosol  size  distribution  and  refractive  index  in  the  sensor  line  of  sight. 

The  atmospheric  constituents  at  any  altitude  will  scatter  the  radiation  incident  on  it;  will 
allow  the  remaining  percentage  of  radiation  to  be  transmitted  onward;  and  also,  will 
emit  a  negligible  degree  of  thermal  radiation  in  this  spectral  region.  The  final  radiance 
that  is  seen  from  a  remote  observer,  looking  either  upwards  or  downwards,  is  the 
accumulated  effect  of  the  transmittance  and  scattering  of  the  atmosphere  on  the  initial 
radiance  {eg.  from  a  target).  The  complications  with  including  the  reflection  off  the 
Earth's  surface  and  accounting  for  this  in  any  retrieval,  suggests  that  downward  looking 
retrieval  via  inversion  techniques  is  unlikely  to  give  accurate  or  stable  results  [19]. 
Particularly  in  the  solar  region,  the  influence  of  the  Earth's  surface  dominates  and 
renders  the  retrieval  of  aerosol  characteristics  almost  impossible. 

Let  us  then  assume  that  we  are  observing  the  atmosphere  along  a  path  that  is  looking 
upwards,  with  a  ground-based  sensor  and  the  primary  radiance  source  from  the  sun.  In 
this  case  we  have  two  scenarios  to  consider,  involving  Sun  and  Sky  radiances. 


4.2  Mathematical  Background 

4.2.1  Aerosol  Size  Distribution  and  Refractive  Index 


The  aerosol  size  distribution  is  a  measure  of  the  number  of  particles  (~  0.001  - 100  pm), 
in  a  unit  volume  of  atmosphere,  which  are  of  a  certain  radii  size  or  fit  within  a  specific 
radii  bin. 


The  function  n{r)  = 


dN{r) 

41ogr 


is  the  aerosol  size  distribution  that  we  wish  to  solve  for.  The 


function  N{r)  is  the  number  concentration  (# / cm^),  which  is  the  number  of  particles  in 
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a  radius  size  bin  per  unit  volume  of  air.  The  total  number  concentration  of  particles, 
then,  is  the  sum  over  all  size  bins,  ^  ^ (^)  •  There  are  four  continuous  size 

distributions  that  are  commonly  used  to  describe  an  aerosol  sample;  the  lognormal,  the 
power  law,  Marshall-Palmer,  and  the  modified  gamma  distributions  [20].  We  will  deal 
mainly  with  the  first  of  these  and  will  thus  henceforth  refer  only  to  the  lognormal  size 
distribution  of  aerosols.  The  appropriate  use  of  other  functions  should  be  considered 
with  reference  to  the  project  application. 


The  lognormal  aerosol  size  distribution  function  is. 


dlogr  ,  V^(log  cr/) 


exp 


-(^logr-logr 

2(logrT;f 


eqn(19) 


where  r  ‘  are  the  mean  aerosol  radii  for  the  population  data,  log  r  ‘  are  the  means  of 
lognormal  distributions,  and  log  <J^  are  the  standard  deviations  of  the  distributions. 

The  complex  refractive  index,  m ,  of  the  atmosphere  is  a  measure  of  the  ability  of  the 
atmospheric  particles  to  absorb  and  scatter  radiation,  m  =  ;/  -I-  .  It  is  a  function  of  the 

wavelength,  m  =  m(A,) ,  and  thus  not  independent  of  the  observed  sensor  radiance 
Aensor  •  The  index  can  be  retrieved  simultaneously  with  the  aerosol  size  distribution  at 
a  fixed  wavelength  [21]  [22],  or,  with  more  difficulty,  as  a  function  of  wavelength  [23] 
[24]. 

As  with  the  retrieval  of  temperature  and  humidity,  the  retrieval  of  both  aerosol  size 
distribution  and  refractive  index  is  achieved  via  reference  to  an  appropriate  radiative 
transfer  equation  (RTE). 


4.2.2  Radiative  Transfer  Equation 

The  scenario  to  be  considered  for  this  application  involves  a  ground-based  electro-optic 
sensor  that  observes  radiation  in  the  solar  region  (0.4  -  l.Sgm)  along  a  path  at  angle  0, 
with  sun  zenith  0^  and  scattering  angles  Q .  With  this  orientation  and  spectral  region, 
we  can  construct  a  RTE  that  neglects  Earth  radiation  and  atmospheric  thermal  emission. 

The  change  in  spectral  radiance  over  the  distance  ds  along  a  path,  gives  the  RTE 

dL(A,  Q)  =  dL^^^^  (A,  Q)  i-  (2-)  —  (A)  —  dL^^^^  (A) 

where  Q  is  the  scattering  angle,  corresponding  to  the  chosen  path  and  sun  orientations. 

Let  us  consider  two  scenarios,  one  in  which  the  sensor  is  directed  at  the  sun  and  one  in 
which  it  is  not.  In  the  first  case,  the  sensor  will  collect  radiation  that  is  coming  directly 
from  the  sun.  The  atmospheric  particles  will  absorb  and  scatter  a  proportion  of  this 
radiation  out  of  the  path,  with  the  remaining  being  transmitted  directly  to  the  sensor.  A 
percentage  of  the  solar  radiation  that  is  not  directly  in  the  line  of  sight  of  the  sensor  will 
also  be  scattered  into  the  path  through  multiple  scattering  processes.  In  the  second  case. 
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when  the  sensor  is  not  directed  at  the  sun,  the  radiance  collected  will  primarily  come 
from  solar  radiation  that  is  scattered  by  atmospheric  particles  into  the  sensor  line  of 
sight  or  'path'.  There  is  no  direct  sun  radiation  component. 

These  two  scenarios  use  the  above  RTE,  with  varying  emphasis  on  the  terms  involved, 
to  gain  information  about  the  atmosphere.  We  will  see  how  this  can  be  combined  to 
retrieve  aerosol  size  distribution  and  refractive  index.  Let  us  call  these  two  cases  direct 
sun  radiation  and  sky  radiation  respectively. 


4.2.2. 1  Direct  Sun  Radiance  RTE 


For  this  scenario  we  have  the  path  zenith  angle  of  the  sensor,  6 ,  equal  to  that  of  the  sun 
zenith  angle,  6^ .  In  this  case,  the  scattering  of  surrounding  radiance  into  the  path  is 
negligible  compared  to  the  direct  solar  radiation  that  will  be  transmitted  (with 
extinction)  directly  to  the  observer.  The  appropriate  approximation  to  the  RTE  includes 
the  effects  of  absorption  and  scattering  extinction  of  the  initial  solar  radiation,  at  various 
wavelengths. 

The  change  in  spectral  radiance  over  the  distance  ds  along  a  beam  is, 
dL{X)  =  (Z)  -  dL^^^^  (A)  +  o{dL^^,^,  (;L))  ;  that  is, 

dL{X)  =  (Z)  -  (Z)  +  o{dL^^^^ ,,,,  (;L))  . 

The  extinction  of  radiation  from  the  sensor's  path  in  the  atmosphere  is  determined  by 
absorption  and  scattering  of  radiation  by  gases  and  aerosols.  The  total  extinction 
coefficient  and  optical  depth  can  be  expressed  as  (A,)  =  (A,)  +  (A,)  and 

^ext  (^)  =  ^sca,  oat  (^)  +  ^absp  (^) '  respectively. 

Therefore,  the  RTE  for  direct  sun  scenario  becomes, 

dL{A)  ~  . 

With  integration  from  space  to  the  ground,  the  direct  sun  RTE  defines  the  radiance 
observed  from  a  ground  based  sensor  at  sO , 

rsO  /  0  \  j  - - 

^sensor  -^RTE  ^o,solar('^)^  -^O.solar 


Expressed  simply. 


rSUN 


=  K 


0, solar 


trans 


atmosphere  ’ 


As  the  initial  solar  irradiance  is  known  for  specific  frequencies  [25],  the  direct 

sun  radiance  is  often  stated  by  the  total  atmospheric  optical  depth,  {A,s0)  =  r‘°[  (A) , 


instead  of  .  That  is; 


(^)  =  Ce  (^,s0)  =  -  cos(0)  log. 


^SUN 

“RTE 


rSUN  /ax 
~^sensor 

V -^0, solar 


eqn(20) 
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With  this  optical  depth  data,  the  direct  sun  radiance  equation  can  be  formed; 


(2)  =  {X,sO)  =  cr  ‘  {X,sO)  +  c 


eqn{21) 


The  optical  depth  r(/l)  is  a  measure  of  the  "penetrability"  of  an  atmosphere  to  a  sensor. 
It  is  an  exponential  function  of  the  transmittance,  which  can  be  used  to  measure  the 
ability  of  the  atmosphere  to  transmit  radiation  that  is  emitted  from  a  remote  source.  The 
various  processes  of  absorption  and  scattering  by  gases  and  aerosols  are  the  main 
influences  that  affect  the  optical  depth  of  an  atmosphere.  Amongst  others,  it  can  be 
shown  that  the  optical  depth  is  a  function  of  the  unknown  variables-  the  aerosol  size 
distribution  and  refractive  index. 


Before  further  defining  the  relationship  between  optical  depth  and  aerosol  size 
distribution  and  refractive  index,  let  us  investigate  the  RTE  of  sky  radiance,  to  draw  a 
comprehensive  picture. 


4.2. 2. 2  Sky  Radiance  RTE 

Consider  the  second  scenario,  in  which  the  path  observation  angle,  0 ,  is  not  the  same  as 
the  sun  zenith  angle,  0^ .  The  upward  looking  observer,  then,  is  not  directly  in  the  path 
of  any  incident  radiation  source,  and  the  only  significantly  contributing  factor  is  the  sun 
radiation  that  is  scattered  into  the  path  and  then  transmitted  (with  extinction)  through 
the  remaining  atmosphere,  to  the  sensor.  The  sky  radiance  RTE  is, 

dL(A,  Q)  =  (A)  +  (A,  Q)  +  {A)) .  eqn{22) 

The  single  scatter  contribution,  ,  is  the  percentage  of  the  solar  radiance  from  the 

sun  that  is  transmitted  to  a  point  on  the  observation  path,  s  * ,  and  then  scattered  into 
the  direction  of  the  sensor; 

‘^Acat  in  =  Aoiar  {A)trans{%w  to  s*)  X  (%  Scattered  into  path)  ds  . 


That  is, 

i.  (-1. n>  =  (l)e  A [o--  (A, Q)  +  C (1)P”" (1, £2)]  * 

where  P  ^  are  the  scattering  Phase  Eunctions  for  gas  or  aerosol  constituents,  as  will 

be  defined  in  the  following  section. 

The  extinction  term,  dL^^^,  in  the  RTE  (eqn(22)),  accounts  for  the  absorption  and 

scattering  of  radiation  out  of  the  remaining  path  to  the  observer,  s*  io  sO .  The  RTE  is 
then  integrated  over  every  scattering  point  along  the  path,  5*  =  5 ,  to  account  for  all 
contributions  to  the  final  observed  radiance. 
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Substituting  the  extinction  and  single  scattering  terms  into  equation  (22),  where 


m(A)  = 


total  /  n  \ 
^scat  (^) 


^  total 


(A) 


is  the  single  scattering  albedo,  the  transformation  of  the  RTE  with 


respect  to  dr  gives; 


dL(A,  Q)  =  L(A,  Q)  +  - 

008(6*) 


^ext 

-^(2)6^"°"^^*^— C7(2) 

°  4;r 


scat 


X=gas,  aero 


gas  +  aero 


P^(A,Q) 


dT(s) 

cos(d) 


eqn(23) 


The  RTE,  in  this  form,  is  not  an  integrable  function  of  optical  depth  or  path  length.  In 
order  to  integrate  the  RTE  and  describe  the  observed  radiance  as  a  function  of  the 
unknown  aerosol  characteristics,  then,  we  must  apply  several  assumptions,  as  detailed 
in  the  following  sections. 


4.2.3  Assuming  a  Homogeneous  Atmosphere 

Consider  an  atmosphere  that  appears,  from  a  ground  sensor,  to  be  optically  equivalent 
to  our  scenario  but  is  in  fact  composed  only  of  a  single  homogeneous  atmospheric  layer, 
of  arbitrary  height.  So  that,  whilst  the  final  optical  depths  are  equivalent  for  this  and  our 
real  scenario,  the  coefficients  of  this  new  scenario  are  now  independent  of  the  altitude. 
By  replacing  our  atmosphere  with  one  such  as  this,  we  are  able  to  integrate  the  sky  RTE 
analytically. 

The  optical  depth  of  the  atmosphere  is  defined  as  the  integral  of  the  coefficient  a(s)  over 

/•oo 

the  path  length;  =  J,  “ cos(^)o-^bsp/ ext/  scat  Ids ' . 

If  we  instead  consider  an  equivalent  atmosphere  that  is  homogeneous  within  an 
imagined  total  path  length,  x  to  5* ,  then  the  corresponding  "uniform"  coefficients,  d", 
are  constant  over  that  path.  The  optical  depth  can  then  be  alternatively  expressed  as, 

Absp/  ext/  scat  (‘^)  “  “  COS(^)  CT^bsp/  ext/  scat  ~  • 

Eor  a  ground  based  sensor,  let  us  define  the  arbitrary  path  length  S ,  such  that 

^absp/  ext/  scat  ^absp/  ext/  scat  * 


Eor  the  sky  radiative  transfer  equation,  this  means  that  the  single  scattering  albedo,  m , 
the  scattering  mixing  ratios  (per  unit  volume)  and  the  aerosol  phase 

function  P’'^°  are  constant  over,  and  thus  independent  of,  the  optical  depth  of  the 
atmosphere. 

With  this  assumption,  we  can  integrate  the  sky  RTE  from  the  top  of  the  atmospheric 
path  to  the  ground  sensor  along  the  line  of  sight,  with  respect  to  optical  depth.  Such  that 
equation  (23)  becomes. 
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f 


cos(0g ) 


cos(0g  )  -  cos(0) 


extW 
COS(^) 


tot 

“^ext 


(A) 


-e 


<x&{6) 


1  _  0 
—  m{A)P  (A,  Q) 
4/r 


J 


where  F^(A,Q) 


<::fw 


P^^\A)  + 


aerosol  /  0  \ 

^scat  )  p  aerosol 

—  total 

T 

scat 


(^,0)  . 


eqn(24) 

eqn(25) 


This  form  of  the  sky  radiance  equation  is  presented  in  terms  of  optical  depth,  rather  than 
coefficients,  since  we  can  express  the  single  scattering  albedo  and  the  scattering  mixing 
ratios  as  follows. 


G7(/l)  = 


^r(^) 

ar(^) 


C‘(^) 


and 


—  gas/  aerosol 
^scat 


(2) 


.  gas/  aerosol 


{X) 


—  total 


{X) 


.  total 


(2) 


eqn(26) 


The  assumption  of  homogeneity  is  in  part  justified  by  the  fact  that  the  spectral  regions 
over  which  we  are  sensing  are  dominated  by  scattering  rather  than  thermal  effects.  The 
ability  of  particles  to  scatter  solar  radiation  is  a  function  of  the  number  and  size  of  the 
particles.  The  majority  of  particles,  especially  larger  aerosols,  are  located  within  the 
lower  altitudes  of  the  atmosphere  (the  boundary  layer).  With  the  orientation  of  our 
sensor  as  ground  based,  looking  upwards,  the  radiance  scattered  into  the  line  of  sight 
will  predominantly  be  the  result  of  the  particles  in  the  lower  altitudes  of  the  atmosphere. 
The  atmosphere  itself  is  obviously  not  vertically  homogeneous,  however,  the  ground- 
based  radiance  can  be  formulated  as  if  it  is,  using  "uniform"  coefficients. 


Thus  we  now  have  the  integrated  sun  and  sky  RTEs  described  as  functions  of  the 
radiance  observed  at  the  sensor  (eqns  (21),  (24)).  Before  we  can  construct  the  relationship 
between  the  observed  radiance  and  the  unknowns,  we  need  to  define  the  various  factors 
involved. 


4.2.4  Assuming  Discretised  Aerosol  Radii  Bins 

Presently,  we  have  sun  and  sky  radiance  equations  that  describe  the  observed  sensor 
data  E^enso,(/i)  as  a  function  of  phase  functions  P^ ,  optical  depth  t^^^(X)  and/or 
extinction  coefficients  (X) .  If  we  wish  to  retrieve  values  for  the  aerosol  size 
distribution  and  refractive  index  from  the  radiance  equations,  we  firstly  need  to  express 
these  variables  in  terms  of  their  constituents.  In  doing  so  we  will  see,  in  the  following 
subsections,  that  these  variables  themselves  are  not  well-posed  and  can  only  be 
described  in  terms  of  discrete  radius  values. 

In  order  to  obtain  well-posed  radiance  functions  for  inversion,  we  need  to  create  discrete 
vectors  to  describe  the  aerosol  size  distribution  and  refractive  index.  Instead  of  a 
continuous  distribution  of  atmospheric  particles,  let  us  construct  logarithmically  spaced 
radii  "bins",  in  order  to  transform  the  aerosol  size  distribution  into  a  vector  format.  The 
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refractive  index  is  a  complex  function  of  wavelength  and  thus  is  represented  by  both  its 
real  ( Re )  and  imaginary  ( Im )  parts,  in  its  vector  form. 

Let  us  define  the  unknown  aerosol  size  distribution  and  refractive  index  as  vector  x , 

X  =  [n^, 

where  n  =  [n(rj),n(r2),...,n(r^)]^and 


4.2. 4.1  Radius-binned  Aerosol  Coefficients  and  Optical  Depths 

The  extinction  and  scattering  coefficients,  ^  ,  measure  the  loss  of  radiation  at  any 

point  along  the  atmospheric  path  due  to  the  presence  of  gases  and  aerosols.  For  aerosols, 
they  are  defined  as  the  product  of  the  respective  cross-sections  and  the  aerosol  size 
distribution.  This  size  distribution  will  vary  through  the  atmosphere,  which  is  often 
represented  by  a  series  of  'sum  of  lognormal'  functions.  Our  assumption  of  vertical 
homogeneity  requires  us  to  treat  this  as  a  vertical  column  of  arbitrary  height  with  a 
single  aerosol  size  distribution  function,  n{r)  . 

C"„,  0.^0}  =  (^>  =  f  = I  (A)  rf(log  r) 

•'  •'  a  log  r 


=  eqn(27) 

where  r  is  the  aerosol  radius,  y^gxt/scat  ^^e  extinction/ scattering  cross-section  for  a 
single  particle,  and  the  function  n(r)  is  the  aerosol  size  distribution  that  we  wish  to 
solve  for. 

The  cross-section  is  defined  as  the  product  of  the  geometric  cross-section  of  a 

single  particle  with  the  extinction/  scattering  efficiency,  ^ext  •  "*"he  efficiency  ^ext  ^ 
complex  function  of  the  wavelength,  particle  radius  and  refractive  index,  m.  The  formula 
for  is  obtained  from  Mie  theory  of  particle  scattering  [26],  [27].  However,  it  can  only 
be  calculated  for  specific  values  of  radii.  That  is,  the  efficiency  is  not  an  integrable 
function  of  radius.  As  such,  the  aerosol  cross-section  y^gxt/scat  integrable  and  the 

aerosol  coefficient  carmot  be  expressed  analytically  as  a  function  of  radius  and  aerosol 
size  distribution.  It  is  necessary,  therefore,  to  discretise  the  aerosol  coefficients  with 
respect  to  radius.  Such  that, 

C7sca,(^n.m)~  z  /^ex,/sca,(^^-."2)-«(^)-A(log^.).  eqn(28) 

i  =radius  bin  1 

The  total  extinction  optical  depth,  is  available  as  observed  spectral  data 

through  manipulation  of  the  direct  sun  RTF,  eqn(20).  With  knowledge  and 


28 


DSTO-TR-1879 


approximations  of  the  molecular  component  of  the  total  optical  depth,  we  can  use  this 
data  to  present  the  direct  sun  total  optical  depth  as  a  function  of  aerosol  size  distribution 
and  refractive  index.  Such  that  the  direct  sun  radiance  equation  becomes. 


i  =radius  bin  1 

eqn{29) 


4.2A.2  Radius-binned  Aerosol  Phase  Function 

The  way  in  which  a  particle  interacts  with  electromagnetic  radiation  depends  on  the 
relative  size  of  the  particle  radius  to  the  radiation  wavelength.  For  gas  species,  the 
wavelength  is  much  larger  than  the  particle  radius  A»  r  and  the  interaction  falls  into 
the  Rayleigh  regime  [28].  This  regime  defines  the  Phase  function  of  such  constituents  to 
be  a  function  of  scattering  angle  alone,  =  fn(Q) . 

For  aerosol  species  however,  the  wavelength  of  solar  radiation  can  be  the  same 
magnitude  as  the  particle  radius  Z  -  r  and  the  interaction  is  described  by  the  Mie 
regime.  In  the  Mie  regime,  the  Phase  function  of  an  optically  homogeneous  atmosphere 
is  defined  as. 


1 

2^ 


j  Amp{A,  Q,  r,  m{X))  n{r)  d  log  r 


where  Amp(A,Q,r,m)  is  a  complex  function  of  the  wavelength,  particle  radius  and 
refractive  index.  The  formula  for  Amp{A,  Q,  r,  m)  can  be  derived  from  the  Mie  theory  of 
particle  scattering  [29],  [30].  Similar  to  the  formation  of  the  efficiency  factor,  however,  the 
Amp  function  can  only  be  calculated  for  specific  values  of  the  radius.  That  is,  it  is  not  an 
analytically  integrable  function,  and  thus  not  in  a  form  suitable  for  analytic  inversion. 

Again,  by  referring  to  radii  "bins",  we  can  discretise  the  phase  function  as  follows; 


A^  1 

^aero(2-,f2,n,m)  =  - - — - - -  ^  Amp{A,Pl,r,,m{A))n{r.)A\ogr.. 

2rr  Acat  n,  jjjj,  [ 


eqn(30) 


With  the  aerosol  phase  function  and  the  aerosol  optical  depth  as  functions  of  n  and  m, 
we  can  now  express  the  sky  radiance  equation  as. 


Fo(A)cos(PJ 
008(6*^  -  008(6*) 


e 


COS(0s) 


V 


cos(^) 


—m(A,m,n)P^(A,Q,m,n) 

4n- 
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where  P®(/l,Q,m,n) 


.  gas 


(Z) 


.  aerosol 


P^^\Z)  +  - 


►  aerosol 


eqn(31) 


We  now  have  equations  for  the  aerosol  coefficients  and  phase  functions  as  discretised 
functions  of  wavelength,  radius,  scattering  angle  and  the  unknown  variables-  aerosol 
size  distribution  and  refractive  index  (eqn(29),  eqn(31)).  With  knowledge  of  the 
molecular  components,  as  given  below,  we  have  constructed  the  observed  sun  and  sky 
radiances  as  functions  of  the  unknown  atmospheric  characteristics.  It  is  these  functions 
that  will  be  the  primary  information  equations  in  the  statistical  regression  and  retrieval. 


4.2.5  Assuming  Molecular  Characteristics 

For  practical  application,  the  observed  radiance  frequencies  are  chosen  in  solar  bands 
that  are  known  not  to  include  molecular  absorption.  Gaseous  absorption  within  the  solar 
spectral  region  is  known  to  be  minimal,  thus  this  is  easily  accomplished; 
^a£!rptioii  ~  0  •  Alternatively,  various  authors  have  accommodated  for  the  presence  of 
ozone  O3  and  nitrogen  dioxide  NO2  in  the  atmosphere  [31],  [32]. 

The  molecular  scattering,  as  described  by  the  Rayleigh  regime,  is  not  dependent  on  the 
size  of  the  gas  species,  as  they  are  more  than  an  order  of  magnitude  smaller  than  the 
solar  wavelengths  that  they  scatter.  An  approximation  of  the  molecular  optical  depth 
can  be  given  via  knowledge  of  the  ground-level  pressure  for  specific  wavelengths  in  the 
solar  region  [33].  According  to  Hansen  and  Travis  (1974),  this  relationship  is, 

Ctout  W  =  0.008569;L  "(1  +  0.01  13;L  "  +  0.000 1  3;l  ")  — 

^0 

where  Z  is  the  wavelength  in  pm,  P  is  the  atmospheric  pressure  in  mb, 
Pq  =1013.25  mb  is  the  standard  pressure. 

4.2.6  "First  guess"  Aerosol  Size  Distribution  and  Refractive  Index 

Similarly  to  the  application  for  temperature  and  humidity  retrieval,  inversions  of  the 
integrated  RTFs  presented  in  equations  (29)  and  (31)  are  not  possible.  Even  though  in 
this  scenario  we  do  have  well-defined  functions  of  the  unknowns,  the  complex 
summation  nature  of  their  structure  makes  them  impossible  to  manipulate.  Thus,  again, 
it  is  recommended  that  we  linearise  the  RTFs  around  an  estimate  of  the  unknowns  and 
construct  the  appropriate  information  equations. 

Fet  us  approximate  the  observed  radiance  about  an  'initial  guess'  profile,  xO,  using  a 
multivariate  Taylor's  expansion. 
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Censor  (^,^)  =  4tE  ~  4tE  ^0)  +  ^ 


M 

z 

k=\ 


5m. 


(m,-m0,) 


xO 


dL^j^{X,Q.) 


5n„ 


+  H.O.T. 


^ 

(n  -nO  ) 

n  n 
xO  J 


In  matrix  form,  the  linearised  equations  are. 


Ls:nL-LRT;(xO)-J^"’'(x-xO)  and  - <te (^0)  -  (x - xO) 


eqn(32) 


where  are  the  Jacobians  of  the  radiance  equations  (eqn(29)  and  eqn(31)). 


jSUN_  ^‘■RTEy 

'*j,i 


5x- 


xO 


^tSKY 

and 

0X, 


eqn(33) 


xO 


and  where  7=1...  {M  +  P)  is  the  number  of  wavelengths  and  scattering  angles,  and 
1  =  1...  {N  +  2M)  is  the  number  of  radii  bins  plus  wavelength  elements. 

Hence  we  have  two  discretised  linear  equations  (eqns(32))  that  require  the  observed 
radiance  and  the  radiance  values  of  the  atmosphere  for  selected  "first  guess" 
atmosphere  conditions  of  n  and  m .  Now  that  the  radiative  transfer  equations  are 
expressed  in  the  generic  form  of  equation  (4)  as  required  in  section  2.2.2,  we  can 
investigate  the  inclusion  of  the  other  information  equations  and  from  thence,  the 
statistical  inversion. 


4.2.7  Retrieval  Equation  for  Aerosol  Size  Distribution  and  Refractive  Index 

For  the  case  in  which  we  are  required  to  retrieve  aerosol  size  distribution  and  refractive 
index,  there  are  four  information  equations  to  consider:  Radiance  and  optical  depth  as 
the  real  data;  and  "first  guess"  and  smoothing  as  the  virtual  information. 

By  combining  the  information  equations  so  that  the  errors  are  minimised,  we  are  able  to 
produce  a  formula  for  the  optimum  solution  to  the  atmospheric  parameter  required,  x  . 
Appendix  A.2.2  has  information  of  the  minimisation  technique  used  and  the  specific 
application  to  aerosol  size  distribution  and  refractive  index  retrieval. 


This  is  now  the  referral  equation  by  which  we  can  calculate  retrieved  values  for 
aerosol  size  distribution  and  refractive  index; 


*»+ (( ism'  is™ ) + (  Jsky'  CEii/' isKY )+ (CeC'  ) + (b'CeC'b))" 


X  = 


T  ^  r 

*'SUN  '-'Elsun 


0 


-1  /  SUN  _-SUN 
^sensor  ^RTE 


(x0)))  +  ( 


T  ^  r 

'^SKY  '^Elsky 


SKY 


-  SKY 


sensor  ^RTE 


(x0)))-( 


B  CEb  BxO 


)) 


eqn(34) 
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where  x  is  the  unknown  aerosol  size  distribution  and  refractive  index;  are 

the  measured  radiances  from  the  sensor;  xO  is  the  approximate  estimate;  Lj^^g(xO) , 
Trte(xO)  are  the  radiances  of  this  estimate;  J  are  the  Jacobians  of  the  radiance 

equations  (29)  and  (31);  and  C  are  the  covariance  matrices  of  the  errors,  as  described  in 
Appendix  B.2. 


4.3  Aerosol  Size  Distribution  and  Refractive  Index  Retrieval: 
Summary  and  Practicalities 

In  summary;  we  consider  the  scenario  in  which  we  are  ground-sensing  solar  radiance  at 
various  wavelengths  and  scattering  angles,  that  has  been  influenced  by  the  intervening 
atmosphere.  We  are  interested  in  retrieving  values  for  the  aerosol  size  distribution  at 
specific  radii  values  and  the  refractive  index  at  specific  wavelengths. 

As  has  been  previously  explained,  the  spectral  data  from  a  remote  sensor  is  not 
sufficient  for  this  problem.  We  must  supplement  the  data  with  a  priori  or  statistical 
information;  for  example,  a  "first  guess"  of  the  distribution  and  index,  from  which  we 
can  extrapolate  the  true  solution.  We  also  require  an  understanding  of  the  errors 
involved  in  the  measurement  and  the  retrieval  procedures,  as  presented  in  the  preceding 
sections. 

Through  statistical  regularisation  of  the  information  available,  the  following  equation 
can  be  obtained,  via  which  we  can  retrieve  the  aerosol  size  distribution  and  refractive 
index  for  a  vector  of  radii  and  wavelengths. 


where. 


X  is  the  vector  of  unknown  aerosol  size  distribution  and  refractive  index  variables  to  be 
retrieved, 

X  =  [n(rj ),...,  n(r^),  Re  (m(,^  )),...,  Re  ( m(4, ))  ,Im(  m(,^  )),...,  Im(  m(4,  ))f ; 


L 

L 


sensor 


sensor 


is  the  vector  of  measured  sky  radiance  values  from  the  sensor, 

“  [Aensor  ’  ^1 )’  Aensor  ^^2  5  ^1 )’  •  •  •  5  ^sensor  ’  ^1 )’  Aensor  (A  ’  ^2  )’•••’  Aensor 


“^sensor  voctor  of  moasurod  direct  sun  optical  depth  values  from  the  sensor. 


t^sensor  Aensor  ('^2  X  *  *  *  5 


xO  is  the  vector  of  first  guess  approximations  of  the  unknown  x , 
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Lj^^g(xO)  is  the  vector  of  sky  radiance  created  by  simulating  an  atmosphere  with  xO 
profile,  using  eqn(31), 

Lrte  (xO)  —  [^RXE  xO), . . . ,  ,  Oj ,  xO),Z,r^e(/^  ,  Qj?  xO), . . .  ,  Qp,  xO)] 

Tr^e  (xO)  is  the  vector  of  direct  sun  optical  depth  created  by  simulating  an  atmosphere 
with  xO  profile,  using  eqn(29), 

"^RTE  (xO)  =  [r 

RTE  (Aj ,  xO), . . . ,  Tr-pe  5  xO)]  } 


J  are  the  Jacobian  matrices  of  the  radiance  equations  (eqn(29)  and  eqn(31)).  The 
derivatives  are  calculated  with  respect  to  the  unknown  variables  and  then  calculated  for 
xO, 


jSUN^  ^*'RTEy 

ax, 


and 


xO 


aySKY 

jSKY_  ^^RTEy 

ax, 


xO 


with  j  =  \...M*  P  frequency  by  scattering  angle  elements  and  i  =  \...N  +  2*  M  radii 
bins  and  frequency  elements; 


and  C  are  the  covariance  matrices  that  account  for  inherent  errors  in  the  process, 

Cei  =  i  [Qi  ,  ^^1 ),  •  •• ,  Qi  and 

Cex  =i[Qx(ni)>-,Qx(nvXC’Ex(miX-,Qx(m2^)f . 


We  initially  assume  that  the  covariance  matrices  are  diagonal.  However,  iteration  during 
the  retrieval  process  can  allow  for  variation,  as  previously  explained  in  section  3.3.I.2. 


4.3.1  Practicalities 

4. 3.1.1  Atmospheric  Models 

Unlike  the  application  of  temperature  and  humidity  retrieval,  the  variables  in  the 
radiative  transfer  equations  that  depend  on  the  unknown  aerosol  characteristics  are  not 
empirically  derived.  That  is,  we  have  calculable  equations  for  the  aerosol  coefficients 
and  phase  function.  Though  these  equations  are  not  continuous  or  simple,  with  the 
discretisation  of  the  RTE  (with  respect  to  radii)  we  have  a  series  of  equations  with 
aerosol  size  distribution  specifically  defined.  The  consequence  of  this  is  that  it  is  possible 
to  find  the  analytic  derivative  of  the  radiance  equation  with  respect  to  aerosol  size 
distribution  n.  This  is  also  true,  but  less  easily  accomplished,  for  the  refractive  index 
m.  Thus  the  Jacobian  matrices  required  in  the  retrieval  equation  above  can  be  obtained 
(Appendices  B.1.2  and  B.1.3). 

The  possible  problems  with  this  approach  include  the  accuracy  of  the  RTE  used  to 
represent  the  atmospheric  processes.  The  radiative  transfer  equations  given  in  the 
sections  above  outline  a  general  approach  to  the  solar  waveband  sensing  scenario,  with 
approximations  made  in  order  to  create  the  retrieval  equation.  However,  more  accurate 
representations  of  the  RTE  can  be  attained  by  using  atmospheric  modelling  software. 
The  choice  to  use  the  analytic  derivatives  to  create  the  Jacobian  matrices,  or  to  use  the 
"small  change"  A  approach  given  in  temperature  and  humidity  application,  is 
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dependent  on  the  specific  circumstances.  For  example;  how  accurately  the  molecular 
contributions  can  be  measured. 

Similarly,  the  first  guess  radiance  (xO)  can  be  calculated  via  forward  modelling  in 
atmospheric  software.  Unlike  the  application  to  temperature  and  humidity  retrieval, 
however,  MODTRAN  does  not  facilitate  the  input  of  new  aerosol  size  distribution  and 
refractive  index  specifically. 


4. 3.1.2  Sources  of  Error  and  Alternatives 

The  retrieval  equation  itself,  in  this  application,  will  need  to  be  dealt  with  carefully  in 
the  calculations.  The  matrices,  especially  the  Jacobians,  if  unchecked,  can  tend  towards 
extremely  large- valued  elements.  The  inversion  of  the  term  required  in  equation  (34) 
then,  which  includes  the  Jacobians,  can  often  lead  to  singularities.  Careful  inclusion  of 

an  appropriate  covariance  matrix  can  reduce  this  possibility.  Various  techniques  for 

mathematically  managing  the  collected  data  and  retrieval  algorithms  are  presented  in 
Dubovik  and  King  [34].  These  include  referring  to  the  radii  in  terms  of  logarithmic 
variables,  and  expressing  the  aerosol  size  distribution  as  a  volume  distribution. 
Supplementing  this,  statistical  techniques  for  confining  the  inversion  term,  in  the 
formation  of  the  relevant  matrices,  can  be  applied  [35]. 

The  retrieval  of  the  refractive  index,  m ,  is  somewhat  more  complicated  than  that  of  size 
distribution.  This  is  due  to  the  fact  that  m(A,)  is  a  complex  function  of  the  coefficients 
and  phase  function,  and  is  dependent  on  the  wavelength  and  thus  not  independent  of 
the  observed  radiance  spectra.  Consequently,  many  applications  of  this  retrieval  will 
treat  the  refractive  index  as  a  constant,  as  this  approach  allows  more  easily  calculated 
Jacobians  and  does  not  require  techniques  to  make  the  solution  distinguishable  from 
noise.  The  paper  by  Romanov  et  al.  [36]  gives  a  novel  approach  to  the  retrieval  problem. 
It  attempts  to  retrieve  refractive  index  as  a  function  of  frequency  and  uses  solar  aureole 
observations  {i.e.  forward  scattered  light  around  the  sun  disk,  which  is  primarily  due  to 
the  scattering  by  aerosols)  to  simplify  the  RTF.  The  use  of  solar  aureole  is  not  unique, 
but  the  retrieval  of  observation  angle  to  make  the  RTE  more  accurate,  is. 

In  regards  to  these  options,  the  chosen  form  in  the  retrieval  of  aerosol  size  distribution 
and  refractive  index  is  dependent  on  the  scenario,  information  available  and  results 
required.  The  retrieval  of  aerosol  characteristics,  from  "solar"  sensors,  can  ultimately  be 
used  to  aid  the  accurate  modelling  of  the  atmosphere  in  a  region,  or  to  aid  the  processes 
of  ground-based  target  detection. 
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5.  Conclusion 


Passive  electro-optic  remote  sensing  can  enable  the  retrieval  of  the  temperature  profile, 
the  humidity  profile,  aerosol  size  distribution  and  the  refractive  index  of  the 
atmosphere.  This  report  has  attempted  to  present  a  generic  design  for  atmospheric 
retrieval,  via  statistical  inversion,  for  the  specific  application  to  sensors  operating  in  the 
long- wave  infrared  and  solar  spectral  wavebands. 

The  inversion  algorithm  is  applied  to  radiative  transfer  equations,  which  have  been 
constructed  using  suitable  mathematical  approximations,  assumptions  and  a  priori 
information.  The  specific  models  of  radiative  transfer  are  given  for  both  of  these 
scenarios  as  a  basic  guide,  but  do  not  encompass  the  complexity  of  the  processes  that  are 
occurring.  With  this  mathematical  foundation,  however,  statistical  analysis  of  the 
retrieval  equations  is  accessible. 

The  inversion  of  the  radiance  equations  is  not  a  well-defined  problem  and  thus  must  be 
solved  via  statistical  techniques  and  regularisation.  Introduction  of  a  "first  guess"  of  the 
unknown  atmospheric  characteristics  redefines  the  problem  as  an  attempt  to  gain  a 
better  solution  to  established  atmospheric  models  or  statistical  a  priori  data.  The 
accuracy  of  the  results,  then,  is  dependent  on  the  error  analysis  and  appropriate 
modelling  of  the  atmosphere  being  considered. 

The  culmination  of  radiance  data  collected  for  both  long-wave  infrared  and  solar  sensors 
could  ultimately  be  used  in  correlation,  to  divulge  information  of  the  gaseous, 
particulate  and  physical  condition  of  the  atmosphere.  This  information,  available  in  the 
temperature  profile,  humidity  profile,  aerosol  size  distribution  and  refractive  index, 
gives  a  more  comprehensive  exposition  of  the  atmosphere.  The  applications  of  this  for 
the  Australian  Defence  Organisation  involve  improved  atmospheric  models  for 
Australian  conditions,  and  more  accurate  atmospheric  compensation  in  target  detection 
algorithms.  To  pursue  these  ends,  the  continuation  of  research  by  DSTO  in  atmospheric 
retrieval  is  recommended. 

This  report  was  written  with  the  aim  to  provide  a  framework  for  atmospheric  retrieval 
via  inversion,  for  a  variety  of  scenarios.  This  framework  can  be  extrapolated  to  include 
phenomena  such  as  clouds  in  the  long-wave  infrared  scenario,  and  applications  such  as 
identification  of  sky-borne  targets.  The  specific  project  and  required  outcomes,  and  the 
availability  of  a  priori  information,  will  determine  the  mode  and  execution  of  the 
atmospheric  retrieval.  Understanding  of  the  assumptions  and  processes  behind  the 
algorithms  used  to  achieve  this  is  vital  for  accurate  variable  representation,  and 
ultimately,  for  reliable  solutions  to  the  atmospheric  characteristics  we  wish  to  retrieve. 

The  purpose  of  this  report  was  the  elucidation  of  the  underlying  mathematics  and 
assumptions  in  the  development  of  a  unified,  generic  description  of  atmospheric 
retrieval.  A  following  report  will  utilise  these  tools  in  application  to  simulated  and 
experimental  data  from  DSTO  scenarios,  with  consideration  for  the  specific  conditions 
present. 
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A.l.  Maximum  Likelihood  Method 

In  Chapter  2,  section  2.3.3,  the  information  equations  are  given  in  the  generic  form, 

4(x)  =  4*+E,, 

where  A:  =  |0,l,2,...if},  4(x)  is  the  true  value  of  the  function  with  the  true  x  value, 
*  is  the  data  that  is  available  or  observed,  and  is  the  uncertainty  or  error. 

The  aim  of  using  the  Maximum  Likelihood  method  is  to  combine  all  data,  both  real  and 
virtual  in  our  scenarios,  to  gain  the  best  estimate  of  the  unknown  vectors.  That  is,  we 
wish  to  maximise  the  conditional  probability  of  finding  f  (x)  given  f  * ; 

P{fl(x),f2(x),...,4(x)  I  . 


There  are  a  few  assumptions  that  need  to  be  made  in  order  to  do  this:  Firstly,  that  all 
information  equations  are  independent  of  each  other;  and  secondly,  that  the  errors  are 
random,  uncorrelated  and  normalised.  If  sufficient  data  is  recorded,  then  the  errors,  if 
unbiased,  will  come  from  a  Normal  distribution  [37].  Using  these  assumptions  enables 
us  to  construct  the  conditional  probability  above  as  a  Gaussian  probability  density 
function; 

P{4(x)  1 4  ~  exp|^(f,(x)-f,  C,'  (f,(x)-f,  *)| 

where  are  the  covariance  matrices  of  the  information  functions. 


And  since,  P{fi(x),f2(x),...,f^(x)  |  *}  =  ]^P{4(x)  |  ^  *},  this  implies 


P{fi(x),f2(x),...,f,,(x)|fi*,f2*,...,f^*}~exp|^^(f,(x)-f,*)^C,'(f,(x)-f,*)l. 


In  order  to  maximise  the  conditional  probability,  we  need  to  maximise  the  exponential 
equation  above.  This  is  accomplished  by  minimising  the  absolute  of  the  index  of  the 
exponential  equation;  that  is,  the  weighted  error  variance. 


A.2.  Retrieval  Equations 

Essentially,  we  need  to  minimise  the  weighted  error  variance  of  the  information 
equations.  In  the  simplest  terms,  the  local  solution  can  be  found  as  follows: 

^  T  - 

Minimising  I(f.  (x)-f^'^)  ^  with  respect  to  x,  implies 

k=0 

-ft,-'  (f,(x)-f,  *)|  =  0  and 


eqn(35) 
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In  order  to  extract  the  unknown  vector  x  from  the  above  equation,  the  particular  form 
of  the  information  equations,  for  the  application  involved,  are  required. 

A.2.1  Application  1:  Temperature  and  Humidity  Retrieval  Equations 
Firstly,  let  us  examine  Application  1;  the  retrieval  of  temperature  and  humidity  profiles, 

x=[ryf. 

Since,  in  equation  (15),  we  converted  the  temperature  and  humidity  radiance  equation 
into  the  general  form  as  described  in  Chapter  2  (eqn  (4)),  we  can  consequently  include 
error  terms  and  present  the  resulting  linearised  equation  as  our  first  information 
equation,  (xO)  +  J  (x -  xO)  ±  El ;  so  that, 

fo*  =  L,ensor  ^ud  fg  (x)  =  (xO)  +  J  (x  -  xO)  , 

where  (xO)  and  J  are  calculated  from  the  temperature  and  humidity  radiance 
equation  (eqn(14)). 


The  second  information  equation  involves  the  error  in  the  "first  guess",  xO  =  [TO^ ,  pO^ 
as  stated  in  the  Chapter  2.3.2,  xO  =  x  ±  Ex ;  so  that, 

fj*  =  xO  and  fj(x)  =  x. 

For  the  retrieval  of  temperature  and  humidity  profiles,  smoothing  is  not  necessary.  Thus 
we  only  need  to  include  these  two  information  equations  in  the  minimisation.  By 
combining  all  information  equations  (x)  such  that  the  errors  are  minimised,  we  are 
able  to  produce  a  formula  for  the  optimum  solution  to  the  atmospheric  parameter 
required  x  in  equation  (35)  above.  Let  us  now  input  these  information  equations  into 
equation  (35)  and  manipulate  to  retrieve  x . 

Firstly,  4-  •}  =  L (xO)  +  j(x  -  xO)  - L„, }  =  j ,  and 

£{f,(x)-f,.}.£{x-xO}.i, 

So  that,  equation  (35)  becomes; 

jd^Er'  +  (CEr')d(L„E(*0>  +  j  (x  -  xO)  -  L.„„, )  +  i  (Ce.-'  +  (Ce.-'  )q(x  -  xO)  =  0 . 

Let  us  assume  here  that  the  errors  are  uncorrelated,  such  that  the  covariance  matrices 
are  diagonal  matrices,  and  thus  C  '  =  (C  *)^ . 

Rearranged,  we  obtain  the  required  retrieval  equation; 
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A.2.2  Application  2:  Aerosol  Size  Distribution  and  Refractive  Index 
Retrieval  Equations 

Let  us  examine  Application  2;  the  retrieval  of  aerosol  size  distribution  and  refractive 
index,  x  =  [n^,m^]^. 


Since,  in  equation  (32),  we  converted  the  radiance  equations  into  the  general  forms 
described  in  the  Chapter  2  (eqn(4)),  we  can  consequently  include  error  terms  and 
present  the  resulting  linearised  equations  as  our  first  information  equations; 


.SUN 


^sVnso.  =  <TE  (xO)  +  j®™  (X  -  xO)  +  El3„^  and  (xO)  +  (x  -  xO)  ±  El 


SKY  ‘ 


So  that, 

=  and  foW  =  T;';;(i:0)  +  j™'‘(x-xO), 
=  and  f,(x)  =  e“(xO)  + j“''(x-xO)  , 


where  Lj^^g(xO)  and  J  are  calculated  from  the  Sun  and  Sky  radiance  equations  (eqn(29) 
and  eqn(31)). 


The  second  information  equation  involves  the  error  in  the  "first  guess"  xO  =  [nO^ ,  mO^ 
which  is  stated  in  Chapter  2.3.2;  xO  =  x  ±  Ex .  So  that, 

fj*  =  xO  and  f2(x)  =  x . 


In  the  case  of  the  aerosol  retrieval,  it  is  also  necessary  that  the  retrieved  solution  be  well 
behaved;  that  is,  smooth  and  continuous.  Let  us  then  introduce  another  virtual 
information  equation  that  will  smooth  the  solution.  As  is  explained  in  Appendix  B.3  the 

appropriate  smoothing  matrix  B  is  chosen  such  that  0  =  Bx  ±  Eb .  So  that, 

f3*  =  0  and  f3(x)  =  Bx. 

For  the  case  in  which  we  are  required  to  retrieve  aerosol  size  distribution  and  refractive 
index  there  are  four  information  equations  ( A:  =  {0, . . . ,  3} )  to  consider:  Sun  radiance  and 
sky  radiance  as  the  observed,  real  data;  "first  guess"  approximation  and  smoothing  as 
the  virtual  information.  Substituting  these  information  equations  into  equation  (35) 
gives; 


( j*” )’  (Ch,„-‘  +  (C,,...-'  f )  (l™  (xO)  +  j*™  ( X  -  xO)  -  )  + 

( + (C„/ V  )(e“  (xO) + j“''  (x  -  xo)  -  , ) + 

i(CHL+(CL)")(x-x0)  +  B"(c,.-‘+(C,.-‘7)Bx  =  0. 

Again,  let  us  assume  here  that  the  errors  are  uncorrelated,  such  that  the  covariance 
matrices  are  diagonal.  That  is,  =  C  so  that  (C)  ’  =  (C^)"'  =  (C  ')^ . 


Rearranged,  we  obtain  the  required  retrieval  equation; 
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Appendix  B:  Retrieval  Matrices 


B.l.  Jacobian  Matrices,  J 

B.1.1  Jabobian  of  "Long-Wave  Infrared"  Radiance 

The  direct  evaluation  of  the  long-wave  infrared  Jacobian  matrix,  corresponding  to  the 
atmospheric  temperature  and  humidity  profiles,  is  complicated  due  to  the  unknown 
composition  of  the  transmittance  function.  An  approximation  to  the  transmittance  can 
be  made  by  substituting  xO  into  an  appropriate  atmospheric  model.  However,  since  we 
do  not  know  the  transmittance  dependence  on  temperature  and  humidity  analytically, 
we  can  therefore  not  calculate  the  required  derivatives  of  transmittance  and  the  radiance 
equation  analytically. 

Theriault  [36]  offers  another  method  for  calculating  J ,  which  treats  the  partial  derivative 
as  a  small  linear  change,  A . 

That  is, 

_  Ly ®  (xO  +  AxO . )  -  (xO) 

AxO, 

xO 

For  the  calculation  of  each  element  of  the  Jacobian  matrix  we  require  the  partial 
derivative  of  the  total  radiance,  with  respect  to  each  layer's  temperature  or 

humidity  value,  evaluated  for  the  "first  guess"  profile  vector.  Theriault's  strategy 
involves  using  FASCODE  or  MODTRAN  to  produce  with  the  original  profile  xO , 
and  then  again  with  profile  xO^,  +  AxO^.  at  layer  k  only.  This  forward-model  run  is 
performed  with  consecutive  changes  to  the  profile,  until  each  layer  has  been  tested 
(Theriault  uses  AxO  =  (2K,0.1g/m^)  [39]).  After  each  "forward  run"  with  an  altered 
layer,  the  total  atmospheric  radiance  is  subtracted  from  the  original  total  radiance,  and 
the  appropriate  Jacobian  is  created  by  dividing  by  AxO  . 

This  approximation  of  the  partial  derivative  matrix  can  potentially  introduce  more  error, 
as  it  assumes  a  degree  of  linearity  of  L’^^Jx)  around  xO.  The  preceding  Taylor's 
approximation  linearised  the  total  radiance  about  xO  for  each  atmospheric  layer,  and 
thus  should  negate  the  effects  of  this  assumption,  depending  on  the  degree  of  the 
linearisation.  The  error  associated  with  linearisation  of  the  RTE  should  be  incorporated 

into  the  covariance  matrix  C^j .  Thus,  AxO^  must  be  sufficiently  small  so  as  to  not  add 
too  large  a  linear  approximation,  but  also  large  enough  so  that  a  sufficient  difference  in 
value  can  be  measured  in  (xO) . 


J,.,.  =  ' 


RTE 


0X, 


AL’ 


RTE 


xO 


Ax,. 


A  combined  numerical  and  analytic  approach  for  calculating  J  is  to  analytically 
differentiate  the  radiance  equations  (eqn(13))  with  respect  to  optical  depth,  and 
numerically  calculate  the  partial  derivatives  of  the  optical  depth  with  respect  to  the 
profiles  at  each  layer.  These  optical  depth  partial  derivatives  would  be  calculated  via 
"  A  "  changes  as  described  above.  This  requires  that  the  optical  depth  can  be  adequately 
represented  by  linear  profile  equations  about  each  altitude,  and  it  this  assumption  that 
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prompts  the  use  of  optical  depth  rather  than  the  transmittance  function  (which  is  an 
exponential  function  of  path  length). 

Again,  the  specific  application,  scenario,  a  priori  information  and  atmospheric  modelling 
packages  available  will  determine  which  method  is  used  to  calculate  the  Jacobian 
matrix. 


^^RTE  ) 

^^RTE  ('^  ) 

^^RTE  (^  ) 

^^RTE  ) 

dT{z^) 

xO 

dT(z^) 

xO 

dp{Zy) 

xO 

dp{z^)  ^ 

^^RTE  i-^M  ) 

^^RTE  (^M  ) 

dd^KYE^^M  ) 

dT{z^) 

xO 

dT{Zj,) 

xO 

dp{Zy) 

xO 

dp{zj,) 

B.1.2  Jacobians  of  "Solar"  Radiance 

The  "solar"  radiance  Jacobians,  required  for  the  aerosol  characteristic  retrieval  equation, 
can  be  calculated  through  analytic  differentiation  of  the  radiance  equations.  This  is 
presented  below.  However,  the  resulting  matrices  can  have  extremely  large  element 
values.  This  can  cause  singularities  during  inversion  in  the  retrieval  equation,  and  thus 
care  in  the  application  of  the  matrices  must  be  taken,  as  specified  in  section  4.3.1. 2. 


For  the  retrieval  of  aerosol  size  distribution  and  refractive  index  we  require  two  Jacobian 
matrices. 


jSUN^  ^^RTEj 


dX; 


and 


jSICY_ 

Jyq  - 


^tSKY 

j 


xO 


dX; 


xO 


For  the  sky  RTF,  is  a  function  of  radius  (bins)  r ,  wavelength  X ,  refractive  index 
m ,  atmospheric  size  distribution  n ,  and  scattering  angle  Q  .  We  require  the  partial 
derivative  of  with  respect  to  size  distribution,  real  part  of  the  refractive  index  and 
imaginary  part  of  the  refractive  index. 


Now,  since, 
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we  can  express  the  integrated  RTE  as, 

(A,Q,sO,m,n)  =  F,  (A)  -  ^ 


1  1 


008(6*5  )  -  cos(6*) 


g  cos(6*s)  cos(6») 


jL71  i  =  bjfj  1 


And  thus,  the  Jacobian  for  the  sky  RTE  with  respect  to  size  distribution,  n ,  is  easily 
calculated  as. 


dL^.^^{A.,D.,sO,m,n)  _ 
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dn{r.)  cos(6*g)-cos(6*) 
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Eor  the  direct  sun  radiance, 

r®TE(^n,m)  =  rf“(;L)+  ^  /3^^^iA,r.,miA)) nir.)  A(logr.) . 

i  =radius  bin  1 


Again,  the  Jacobian  elements  are  easily  calculated  with  respect  to  the  aerosol  size 
distribution. 


arRT^(;L.,5d,n,m) 

dn{r.) 


=  /^ext  (A  ’  n  ’  ))  ^(log  n )  • 


To  calculate  the  analytic  derivatives  of  the  radiance  equations  with  respect  to  aerosol 
size  distribution,  therefore,  we  required  the  Mie  formulae  for  Amp(A,r)  and  /3(A,r)  . 
These  complex  summation  formulae  are  specified  in  Deirmendjian  [39]  for  example. 

Eor  the  Jacobians  of  the  refractive  index,  the  derivatives  are  not  as  simple  to  determine 
analytically.  That  is,  since  Amp  is  a  function  of  m ,  the  partial  derivatives  of  the  real  and 
imaginary  parts  of  the  refractive  index  are  non-trivial. 

Eor  the  real  part  of  the  complex  refractive  index  function,  the  Jacobian  elements  are 
described  by, 
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Similar  matrix  elements  can  be  derived  for  the  imaginary  part,  Im(m) ,  for  both  sun  and 
sky  Jacobians. 

A  useful  construction  of  the  derivatives  of  the  Amp{X,  r)  and  /3{A,  r)  functions  with 
respect  to  the  refractive  index  is  given  in  [41]  and,  with  greater  detail,  in  [42]. 


Below  is  the  sky  Jacobian, 
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A  similar  Jacobian  matrix  can  be  obtained  for  the  direct  sun  radiance  (without  reference 
to  the  scattering  angle,  Q ). 


B. 2.  Covariance  Matrices,  C 

The  covariance  matrices,  ,  are  introduced  as  the  error  variance  parameters  in  the 
probability  formulae  used  to  optimise  the  solution  (Appendix  A.l).  The  diagonal 
elements  of  are  defined  by  the  square  of  the  error  in  the  information  equations,  . 
It  is  through  setting  the  covariance  parameters  suitably  that  the  solution  can  be 
constrained  (section  3.3.1). 

The  application  of  probability  theory  is  used  to  maximise  the  probability  of  obtaining 
f^(x)  given  *,  in  order  to  extract  an  accurate  value  for  x  with  minimum  error.  We 
define  this  error  as  ,  such  that  4  *-  f^-Jx)  =  E^ .  The  error  is  simply  the  difference 

between  the  true  solution  and  the  approximation.  The  diagonal  elements  of  the 
covariance  matrices  are  then  defined  by  the  appropriate  error  variances, 

C, (/,0  =  E,(0-E,(0. 

The  covariances  are  firstly  assumed  to  be  the  result  of  independent  errors,  such  that  the 
matrices  are  diagonal  {i.e.  all  non-diagonal  elements  are  zero).  With  Newtonian  iteration 

of  the  solution,  it  is  also  possible  to  iteratively  improve  the  covariance  matrix,  ,  for  a 

more  accurate  solution.  See  [43],  [44]  and  [45]  for  greater  mathematical  detail,  in 
retrieving  different  scenarios. 
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For  the  first,  real  information  equation  (xO)  +  J(x  -  xO)  ±  El ,  so  that 

the  variance  is  (Lsensor-(LRTE(xO)  +  J(x-xO))j  =  (±El)" .  The  diagonal  covariance 

matrix  is  then  defined  by  Cgj (/,  i)  =  El(/)  •  El(/) .  As  can  be  similarly  shown  for  total 
optical  depth,  . 

The  error  El  is  a  combination  of  the  sensor  error  and  the  error  in  the  linearised  radiance 
equation,  El  =  .  These  values  must  be  estimated  from  the  known  sensor 

sensitivity  and  variability,  and  the  consideration  of  the  a  priori  model  known  for  the 
atmosphere.  Setting  an  appropriate  covariance  matrix  Cg[  can  help  stabilize  the  retrieval 
by  including  the  "actual  limit"  of  the  measurement  and  the  uncertainty  in  the 
atmospheric  model. 

For  "first  guess"  virtual  information  we  have  fj  *  =  xO  =  x  ±  Ex ,  so  that  the  variance  is 
(xO  —  x)^  =  (±Ex)^ .  The  diagonal  covariance  matrix  is  then  defined  by 
CEx(hO  =  Ex(0-Ex(0. 

The  covariance  matrix  of  the  first  guess  of  the  solution,  ,  is  often  referred  to  as  the 
"damping  factor".  The  elements  of  the  matrix  are  arbitrarily  chosen  in  order  to  stabilise 
the  retrieval  by  confining  the  solution  to  within  a  chosen  "error"  degree  of  the  "first 
guess".  A  better  estimation  of  the  covariance  elements  can  then  be  gained  by  iterating 
the  solution,  as  described  above.  This  is  not  possible  for  all  scenarios,  and  since  the 
covariance  matrix  is  ultimately  determined  by  "chosen"  values,  this  parameter  of  the 
retrieval  equation  is  not  a  robust  variable. 

For  the  smoothing  we  have,  f4*  =  0  =  Bx±Eb,  so  that  the  variance  is 
(0-Bx)"  =(±Eb)"  .  The  diagonal  covariance  matrix  is  Cgg(/,/)  =  Eb(0  •Eb(z) . 


B.3.  Smoothing  Matrix,  B 

The  retrieved  solution  for  the  aerosol  size  distribution  is  known  to  be  a  smooth 
function.  This  'smoothness'  can  be  used  as  virtual  information  to  obtain  a  more  accurate 
inversion  [46].  The  degree  to  which  the  function  is  to  be  smoothed  determines  the  form 

of  B .  The  curvature  of  a  function  is  defined  by  its  second  derivative,  and  thus 
restrictions  on  this  can  influence  the  smoothness.  That  is,  by  minimising  the  curvature  or 
the  change  in  curvature  (third  derivative),  the  function  should  appear  to  become 

smoother.  Dubovick  and  King  [47]  describe  the  smoothing  parameter,  B ,  as  a  matrix 
operator  of  "second  differences"  on  the  profile  x  [ie.  the  coefficients  of  recurrent 
difference  equations,  as  an  approximation  to  derivative  functions).  Within  the  inversion 

equation  it  is  this  product,  B  x ,  that  is  minimised.  That  is,  we  require  the  second 
derivative  of  the  solution  to  tend  to  zero,  with  minimum  error,  and  thus  give  a  smoother 
solution,  X . 
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